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Summary  of  Completed  Project 

This  AASERT  program  supported  research  in  robust  simulation,  hierarchical  uncertainty 
representation,  and  novel  methods  for  robustness  analysis  of  uncertain  systems. 

In  the  context  of  this  program,  robust  simulation  means  simulating  simultaneously  sets  of 
initial  conditions  and  disturbance  or  noise  signals.  Thus  sets  of  state  space  must  be  propogated 
by  the  dynamics  of  the  model.  Initial  investigations  have  focused  on  piecewise  linear  discrete 
time  systems,  which  map  polyhedra  to  polyhedra  at  each  time  step.  Linear  programming  can 
be  used  to  approximation  complicated  polyhedra  with  simpler  bounds,  and  branch  and  bound 
can  be  used  to  refine  the  resulting  bounds.  This  is  important  if  the  potentially  exponential 
growth  in  set  descriptions  is  to  be  overcome. 

Hierarchical  uncertainty  modeling  is  a  new  framework  to  include  explicit  representation 
of  uncertainty  in  component  modeling.  The  focus  has  been  on  LFTs  and  implicit  (DAE) 
representations.  A  variety  of  examples  including  parasitics  and  nonlinearities  illustrate  the 
key  ideas. 

Finally,  this  report  describes  new  bounds  on  a  “spherical  /j, ”  problem  that  allows  for 
correlations  between  uncertainties  in  an  LFT  framework.  Interestingly,  this  setting  provides 
quite  elegant  bounds  and  simplified  computation. 
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1  Robust  Simulation 

1.1  Introduction 

Historically,  simulation  has  played  a  large  role  in  nonlinear  system  analysis.  Stability  and 
performance  are  often  studied  by  simulating  a  large  number  of  initial  conditions  and  noise 
signals.  However,  the  results  of  these  simulations  do  not  guarantee  stability  or  performance, 
since  the  behavior  for  other  initial  conditions  and  noise  signals  is  not  examined. 

Robust  simulation  addresses  this  limitation.  Instead  of  simulating  a  single  initial  condition 
and  noise  signal,  a  set  initial  conditions  and  noise  signals  are  simulated  at  once.  If  the 
robust  simulation  meets  the  performance  requirement,  then  the  system  meets  the  performance 
requirement. 

As  with  any  general  nonlinear  problem,  some  assumptions  are  required  to  ensure  reason¬ 
able  computational  cost.  For  robust  simulation,  finite  time  horizon  problems  for  discrete  time 
piecewise  linear  systems  are  considered.  All  simulation  based  techniques  require  finite  time 
horizons  and  the  discrete  time  piecewise  linear  restriction  is  less  onerous  than  it  sounds. 

Traditional  simulation  maps  a  single  point  in  initial  condition  and  noise  space  into  a 
single  final  condition.  We  define  robust  simulation  as  the  mapping  of  all  initial  conditions 
and  noise  signals  into  a  set  of  final  conditions.  By  calculating  all  possible  trajectories,  one  can 
make  guarantees  about  the  system’s  global  behavior.  Traditional  simulation  only  gives  local 
information.  Essentially,  robust  simulation  maps  sets  to  sets,  while  traditional  simulation 
maps  points  to  points. 

The  mapping  of  sets  introduces  a  new  problem,  set  representation,  into  the  simulation 
process.  For  linear  systems,  this  is  not  a  difficult  issue,  since  polyhedra  map  to  polyhedra 
via  matrix  multiplication.  However,  nonlinear  systems  can  map  polyhedra  into  anything.  For 
efficient  computation,  simple  set  descriptions  and  mapping  rules  are  required. 

The  set  representation  issue  is  simplified  by  requiring  discrete  time  piecewise  linear  sys¬ 
tems,  which  map  polyhedra  to  polyhedra  via  simple  matrix  operations.  This  simplification 
also  eliminates  two  problems  encountered  in  traditional  simulation:  step  size  determination 
and  derivative  approximation.  This  allows  us  to  focus  on  the  issues  central  to  robust  simula¬ 
tion,  set  representation  and  computational  complexity. 

In  order  to  simulate  sets,  assumptions  must  be  made  on  the  form  of  system.  By  considering 
only  PL  systems,  all  mapping  operations  can  be  accomplished  through  linear  programming 
and  simple  matrix  operations.  Since  large  linear  programs  are  easily  solved  by  existing  soft¬ 
ware  and  linear  program  size  grows  reasonably  with  state  dimension,  large  systems  can  be 
reliably  simulated. 

The  idea  of  the  algorithm  is  quite  simple.  Start  with  an  initial  condition  set  and  noise 
description  and  calculate  all  possible  final  conditions  after  one  time  step.  This  is  repeated 
until  the  desired  final  time  is  reached.  While  this  simple  approach  calculates  the  exact  set 
that  can  be  reached  at  the  final  time,  it  also  has  exponential  computational  cost.  Assuming 
that  the  initial  set,  <So,  is  a  single  polyhedra,  the  reachable  set  at  any  time  £,  denoted  S*,  may 
contain  lf  polyhedra.  Since  So  may  intersect  with  each  of  the  it!*,  different  subsets  of  <Sq  may 
follow  up  to  l  different  state  update  laws,  each  yielding  a  separate  polyhedra.  Thus,  <Si  may 
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contain  l  polyhedra.  Each  polyhedra  in  S\  may  map  into  /  additional  polyhedra,  yielding  up 
to  / 2  polyhedra  in  £2-  For  long  simulations,  this  exponential  growth  is  unacceptable. 

By  restricting  the  number  of  polyhedra  in  each  set  Sk,  exponential  growth  is  avoided.  This 
is  accomplished  by  combining  polyhedra  into  larger  polyhedra  such  that  each  of  the  original 
polyhedra  is  a  subset  of  the  new  one.  Thus,  Sk  contains  all  points  reachable  from  Sk-i  plus 
additional  points  added  when  combining  polyhedra.  This  approximation  is  conservative.  The 
simulation  result  always  contains  all  reachable  points,  but  may  contain  many  points  that  are 
nor  reachable.  However,  the  conservatism  can  be  systematically  reduced  by  allowing  more 
polyhedra  in  Sk-  When  Sk  can  contain  lk  polyhedra,  the  exact  answer  is  obtained. 

With  the  above  approximation,  robust  simulation’s  computational  cost  grows  linearly  with 
time.  At  each  time  step,  the  same  number  of  polyhedra  are  being  mapped.  Doubling  the 
simulation  length  only  doubles  the  computation  time.  Computational  cost  grows  polynomially 
in  other  variables,  as  well.  The  approximation  requires  solving  roughly  tnP+2  linear  programs, 
where  0  <  7  <  t  is  the  degree  of  refinement.  For  basic  robust  simulation  (7  =  0)  using  a 
public  domain  0(ri 4)  linear  program  solver,  the  overall  complexity  is  0(tl2n 5).  Memory  usage 
is  0{llJt2nd).  These  are  worst  case  complexities;  performance  is  generally  much  better  than 
0{tl2n 5). 

1.2  Piecewise  Linear  Systems 

When  developing  numerical  tools  for  nonlinear  analysis,  piecewise  linear  (PL)  systems  are 
a  natural  class  of  systems  to  consider.  They  are  a  conceptually  simple  extension  of  linear 
systems.  In  any  region  of  state  space,  an  affine  state  update  law  is  followed.  PL  systems  are 
well  suited  to  digital  implementation  and  simulation;  only  matrix  manipulations  and  if-then 
blocks  are  needed.  PL  systems  can  be  completely  described  by  a  finite  amount  of  information. 
Even  the  most  complex  PL  system  is  described  by  a  list  of  matrices. 

As  simple  as  they  sound,  PL  systems  are  a  very  rich  class.  Many  nonlinearities,  such  as 
saturations,  dead  zones,  quantizations,  and  hysteresis,  are  naturally  written  as  PL  mappings. 
PL  systems  can  exhibit  very  complex  dynamics,  including  chaos.  In  fact,  PL  systems  can 
approximate  general  nonlinear  systems  to  any  degree  of  accuracy. 

Unfortunately,  the  simple  intuition  and  ease  of  implementation  does  not  extend  to  analysis. 
Calculating  almost  any  interesting  PL  system  property  is  NP-hard.[35]  Decades  of  study 
have  yielded  few  practical  analysis  results.  In  1981,  Sontag  suggested  the  use  of  PL  systems 
for  nonlinear  regulation  [33]  and  developed  a  PL  algebra  [34].  Recently,  Pettit  and  others 
have  studied  continuous  time  PL  systems. [31]  Unfortunately,  these  results  have  not  yielded 
practical  analysis  tools.  Computations  either  grow  exponentially  with  problem  size  or  do  not 
give  guarantees. 

A  piecewise  linear  (PL)  system  is  defined  over  some  subset,  Z,  of  a  finite  dimensional  real 
vector  space,  1Zn.  Z  is  the  union  of  a  finite  number,  l,  of  closed  polyhedra,  denoted  Ri .  Each 
Ri  is  defined  by  a  finite  number  of  linear  inequalities,  f(x)  >  a.  In  each  i?;,  an  affine  state 
transition  map  is  defined  by 

x(k  +  1)  —  Aix(k)  +  Bi  +  Niw[k\,  x  e  Ri,  i€l...i.  (1) 
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The  noise  signal,  w[k\,  is  chosen  to  be  in  to  simplify  computations.  The  methods  presented 
trivially  extend  to  time-varying  PL  systems,  as  well. 

This  definition  simplifies  computer  implementation,  but  leads  to  a  minor  technical  prob¬ 
lem.  Since  the  polyhedra  are  closed  and  may  share  boundaries,  the  mapping  is  not  always 
well  defined  on  the  boundaries.  While  it  is  possible  to  create  a  PL  system  that  exploits  the 
behavior  on  the  boundaries,  this  is  generally  not  the  case.  Because  this  technical  issue  is 
easily  resolved  (see  [33]),  and  only  clutters  the  description  of  the  problem,  it  will  be  ignored 
for  the  remainder  of  this  presentation. 

By  construction,  each  Ri  is  convex  and  bounded  by  hyperplanes.  While  the  number  of 
hyperplanes,  c,  bounding  each  Ri  may  vary,  c  must  be  greater  than  n  +  1  (to  define  a  simplex) 
and  generally  c  —  2n  (to  define  a  hyperrectangle).  The  notation  Sj  denotes  a  finite  set  of 
polyhedra  in  TZn  at  time  j. 

Linear  programming  and  convex  set  theory  play  a  large  role  in  PL  system  analysis.  By 
construction,  each  volume  Ri  is  convex  and  bounded  by  hyperplanes.  While  the  number  of 
hyperplanes,  c,  bounding  each  region  Ri  may  vary,  c  must  be  greater  than  n  +  1  (to  define  a 
simplex)  and  generally  c  =  2n  (to  define  a  hyperrectangle).  The  hyperplanes  bounding 
called  Hi ^  U1...C,  are  defined  by  a  point,  Xifa  and  a  normal  vector,  Nifa  The  normal 
vector  is  defined  as  pointing  into  the  region  Ri.  With  this  definition,  linear  programming  can 
be  used  to  solve  a  variety  of  problems.  For  example,  a  point  x  lies  inside  Ri  if  and  only  if 

(x  -  xith  ■  Nij l)  >  0,  V/i  G  1 . . .  c 

This  can  be  written  as  a  linear  program  with  c  constraints.  Other  problems,  such  as  finding 
the  intersection  of  regions  and  determining  norms  on  regions  are  also  linear  programs. 

1.3  Robust  Simulation  Algorithm 

The  robust  simulation  algorithm  answers  the  following: 

For  a  given  set  of  initial  conditions,  denoted  So,  what  are  all  possible  final  condi¬ 
tions,  denoted  Sy? 

By  assigning  a  norm  on  Sy,  this  becomes  a  performance  measure. 

The  direct  approach  to  calculating  Sy  has  exponential  growth  in  the  computation  as  a 
function  of  t.  To  see  this,  start  with  So,  map  forward  one  time  step,  and  call  the  new  set  Si. 
Assuming  (for  notational  simplicity)  that  So  contains  at  most  Z  convex  regions,  there  are  at 
most  Z 2  convex  sets  in  the  So  CiRi.  Since  no  restrictions  are  placed  on  (1),  each  of  these  convex 
sets  can  then  map  into  all  of  the  Ri.  Thus,  Si  can  contain  up  to  Z2  convex  sets  and  Si  D  Ri 
may  contain  Z3  sets.  At  any  time  step  j,  Sj  can  contain  as  many  as  ZJ+1  independent  convex 
sets.  Repeating  this  process  to  form  St  yields,  possibly,  Zt+1  sets.  All  of  these  mappings  can 
be  computed  by  simple  matrix  operations  and  linear  programming. 

By  slightly  modifying  the  direct  approach,  a  polynomial  time  bound  is  obtained.  The 
fundamental  idea  is  to  limit  the  number  of  regions  in  Sj  at  every  time  step  by  restricting 
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Sj  to  have  a  fixed  number  of  convex  sets  in  each  R{.  The  restrictions  placed  on  Sj  deter¬ 
mine  the  tightness  of  the  bound.  The  result  contains  Sf,  though  the  approximation  may  be 
conservative. 

Given  Sj,  the  first  step  is  to  form  a  manageable  approximation  Tj  2  Sj.  <Sj+i  is  then 
calculated  from  Tj.  By  restricting  the  number  of  sets  in  the  approximation,  exponential 
growth  is  avoided.  By  definition,  Tj  contains  P+1  convex  sets,  and  there  are  V  convex  sets 
in  Tj  n  Rj,.  The  meaning  of  7  will  be  described  later.  Though  the  details  are  omitted  due  to 
space  constraints,  Tj  can  be  calculated  by  linear  programming.  [12] 

1.3.1  Algorithm  Refinements 

The  accuracy  of  the  approximation  is  directly  related  to  the  amount  of  extra  volume  added 
when  forming  Tj.  Two  factors  affect  this:  the  number  of  convex  sets  in  each  Ri,  and  the 
number  of  hyperplanes,  c,  used  to  bound  each  region.  By  increasing  each  of  these,  accuracy 
may  be  improved. 

As  defined  earlier,  Tj  fl  R,  contains  P  convex  sets.  The  value  7  can  be  considered  a  history 
parameter.  For  7  =  0,  one  does  not  consider  what  regions  a  set  mapped  from  before  arriving 
in  the  current  Ri.  For  7  =  1,  one  looks  at  where  the  set  was  one  time  step  prior  to  the  current 
time  step.  In  this  case,  each  Tj  n  Ri  contains  l  convex  sets,  each  approximating  the  sets  that 
came  from  a  specified  R.  7  determines  how  many  previous  time  steps  play  a  role  in  forming 
the  approximation  Tj.  As  7  approaches  t,  the  approximation  approaches  the  exact  solution. 
When  7  =  t,  the  exact  solution  is  obtained. 

An  equally  important  variable  is  the  number  of  hyperplanes  defining  a  convex  set,  c. 
Assuming  hyperrectangles,  the  basic  algorithms  required  2 n  bounding  surfaces.  In  general, 
this  is  2n(7  + 1).  The  bounding  hyperplanes  must  contain  those  of  the  current  region,  Ri,  and 
the  bounds  of  the  previous  regions  considered  after  the  mapping  law  is  applied.  This  way,  a 
hyperrectangle  can  be  exactly  covered  after  mapping  7  time  steps. 

Since  c  >  2n(7  +  1),  an  ad  hoc  approach  for  improving  the  approximation  is  to  add  more 
bounding  surfaces  according  to  some  heuristic.  In  general,  additional  hyperplanes  should 
differ  greatly  from  existing  hyperplanes. 

Additionally,  other  generic  bound  improvement  techniques,  such  as  branch  and  bound  are 
applicable. 

1.3.2  Computational  Complexity 

The  robust  simulation  approximation  requires  solving  tcP+ 2  linear  programs.  0(n4)  public 
domain  linear  program  solvers  are  available  and  0(nz-b)  algorithms  have  been  proposed.  [17] 
With  c  =  2 n,  the  default  when  using  hyperrectangles,  7  =  0,  and  an  C>(n4)  linear  program 
solver,  the  overall  complexity  is  0(tl2n5). 

Memory  usage  is  0(P+2nc).  At  any  time,  up  to  P+2  linear  programs  must  be  stored  in 
memory.  These  linear  programs  require  memory  proportional  to  nc.  Since  the  results  from 
previous  time  steps  do  need  to  be  saved,  memory  usage  is  independent  of  t. 
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Figure  1:  Nonlinear  performance  bound  ratios 

In  general,  performance  is  much  better  than  0(tl2nb).  This  worst  case  performance  as- 
sumes  that  each  region  Ri  maps  into  every  other  region  at  each  time  step.  Many  systems 
map  only  into  adjacent  regions  or  themselves  at  each  time  step.  By  calculating  what  regions 
flow  into  other  regions  in  advance,  the  number  of  linear  programs  solved  at  each  step  can  be 
greatly  reduced. 

1.3.3  Examples 

To  evaluate  the  algorithm,  tests  were  run  on  a  set  of  randomly  generated  systems.  The  random 
systems  were  discretizations  of  5th  order  continuous  systems  with  one  saturation  nonlinearity 
and  no  noise.  Systems  were  simulated  for  30  time  steps.  The  performance  measure  chosen 
was  | x  [30]  |oo.  The  initial  condition  set  was  |^[0]|oo  <  1.3.  This  is  a  reasonable  set  of  test 
problems  that  are  moderately  challenging  for  the  algorithm.  This  is  neither  the  hardest  nor 
the  easiest  class  of  problems  known. 

Any  traditional  simulation  gives  a  lower  bound  on  the  worst  case  performance.  To  achieve 
a  large  lower  bound,  gradient  search  was  combined  with  random  simulation.  The  upper  bound 
was  calculated  by  robust  simulation  with  7  =  0.  If  the  bounds  differed  by  more  than  10%, 
naive  branch  and  bound  was  applied  until  the  bounds  differed  by  less  than  10%  or  50  branch 
steps  were  taken. 

The  measure  of  algorithm  performance,  shown  in  Figure  1,  is  the  ratio  of  the  lower  bound 
to  the  upper  bound.  Ideally,  the  ratio  is  always  1.  In  general,  ratios  greater  than  0.9  are 
acceptable.  For  89%  of  the  runs,  no  branching  was  needed  to  obtain  acceptable  results.  For 
the  11%  of  the  runs  where  branching  was  needed,  branch  and  bound  greatly  improved  the 
bound  ratio.  These  results  can  be  greatly  improved  by  better  selection  of  branch  cuts. 

1.4  Robust  Stability  of  PL  Systems 

Efficient  algorithms  for  verifying  the  stability  of  uncertain  discrete  time  piecewise  linear  (PL) 
systems  will  be  presented.  While  PL  systems  are  intuitively  simple,  they  are  computationally 
hard.  Two  approaches  to  verifying  stability  are  presented.  For  each  approach,  separate 
necessary  and  sufficient  conditions  are  given. 
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The  first  approach  requires  the  solution  of  a  linear  matrix  inequality.  This  method  is 
only  applicable  to  a  restricted  class  of  PL  systems,  and  is  generally  very  conservative.  It  is 
demonstrated  that  for  most  PL  systems,  these  conditions  yield  no  information. 

The  second,  more  general,  approach  is  based  upon  robust  simulation.  This  method  is 
useful  for  all  PL  systems,  and  will  always  yield  a  definitive  answer.  If  a  system  initially 
satisfies  necessity,  but  fails  sufficiency,  these  algorithms  can  be  systematically  refined  and 
after  a  finite  number  of  refinements,  a  definitive  answer  is  guaranteed. 

The  algorithms  are  illustrated  on  four  examples,  including  examples  where  no  other  general 
analysis  technique  is  available. 

1.4.1  LMI  based  Alogorithm 

For  a  restricted  class  of  PL  systems,  separate  necessary  and  sufficient  conditions  for  stability 
require  solving  LMI’s.  However,  these  conditions  are  very  conservative. 

For  PL  systems  without  any  affine  terms  ( Bi  =  0  Vi),  stability  can  be  guaranteed  by 
finding  a  feasible  solution,  V,  to  the  LMI 

V  >  0 

A\VAi  —  V  <  0,  *  €  1 ...  I.  (2) 

Similarly,  instability  is  guaranteed  if  a  feasible  solution  to 

V  >  0 

A\V Ai  -  V  >  0,  i  G  1 . . .  I  (3) 


exists. 

If  (2)  has  a  solution,  then  a; [A;]* Fa; [A:]  is  a  common  quadratic  Lyapunov  function  for  all 
Ai  and  the  PL  system  is  stable.  This  Lyapunov  function,  however,  is  independent  of  the 
switching  law.  Switches  can  occur  arbitrarily  quickly  or  even  randomly,  and  the  PL  system 
is  still  stable.  Clearly,  this  is  sufficient  condition  for  stability,  though  it  is  not  necessary.  In 
section  1.4.3,  an  example  is  given  that  fails  this  test,  but  is  still  stable. 

For  stability,  it  is  also  necessary  that  there  is  no  solution  to  (3).  If  a  solution  exists,  then 
V(A;)  =  a; [A;]* Fa; [A;]  is  Lyapunov-like  function  where  V(A;  +  1)  >  V(A;)  VA;,  and  the  system  is 
unstable  for  any  switching  law. 

As  previously  mentioned,  these  conditions  are  very  conservative.  If  any  Ai  has  an  unstable 
pole,  then  (2)  never  has  a  feasible  solution.  One  can  easily  construct  a  stable  PL  system  where 
some  A{  are  unstable.  If  any  A{  is  stable,  then  (3)  never  has  a  solution.  Section  1.4.3  presents 
a  PL  system  where  all  Ai  are  stable,  but  the  system  is  unstable.  These  results  are  closely 
related  to  discrete  linear  inclusions  and  are  well  known.  [3,  10] 

1.4.2  Robust  Simulation  Algorithm 

The  robust  simulation  based  algorithm  examines  generalized  finite  time  horizon  stability, 
specifically: 
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Do  all  states  in  a  set  of  convex  regions  Sq  map  into  a  set  of  convex  regions  Sf 
within  t  time  time  steps? 

With  a  careful  choice  of  Sf,  stronger  notions  of  stability  can  be  inferred.  If  Sf  is  known 
to  be  invariant,  then  Lyapunov  stability  is  implied.  If  Sf  is  known  to  be  invariant  and  a 
quadratic  Lyapunov  function  exists  for  Sf,  then  quadratic  stability  is  implied.  The  latter 
case  occurs  frequently,  often  when  0  C  Sf  C  Ri  for  a  fixed  i.  Answering  this  question  is  a 
direct  application  of  robust  simulation,  which  calculates  all  possible  trajectories  for  a  set  of 
initial  conditions. 

A  sufficient  condition  for  generalized  stability  is  obtained  by  robustly  simulating  So  for 
t  time  steps,  yielding  the  set  of  all  possible  final  points,  Sf.  If  St  C  Sf,  then  the  stability 
condition  is  satisfied.  It  should  be  noted  that  robust  simulation  is  conservative;  St  contains  all 
possible  final  conditions  and  additional  points  that  cannot  be  reached  from  So.  This  condition 
only  guarantees  stability;  if  it  is  not  met,  nothing  can  be  inferred. 

A  necessary  condition  for  generalized  stability  is  found  by  robustly  simulating  Sf  for  t 
time  steps  backwards  in  time,  yielding  the  set  Sf-t-  Robust  simulation  in  backwards  time 
is  always  possible,  even  if  some  A{  are  not  invertible.  If  So  %  Sf-t ,  then  the  system  not 
stable.  When  Sq  %  Sf-t,  some  point  in  <So  does  not  map  into  Sf  after  t  time  steps.  Since 
Sf-t  contains  points  that  do  not  map  into  Sf ,  this  is  only  a  necessary  condition. 

If  the  sufficient  condition  fails  and  the  necessary  condition  is  satisfied,  no  conclusion  can 
immediately  be  drawn.  However,  the  conditions  can  be  refined  by  using  more  accurate  robust 
simulation  algorithms  at  the  expense  of  computational  cost.  The  exact  answer  can  always  be 
obtained  in  finite  time. 

The  generalized  stability  problem  can  easily  be  rephrased  as  an  3  problem  by  checking 
for  instability.  For  PL  systems,  3  problems  are  in  general  NP-hard.[35]  Thus,  an  efficient 
algorithm  that  gives  the  exact  answer  is  not  expected.  The  question  can  be  shown  false  in 
polynomial  time  by  choosing  an  appropriate  xq  and  calculating  is  trajectory. 

These  algorithms  extend  directly  to  uncertain  (noisy)  PL  systems.  Robust  simulation  can 
be  used  to  analyze  systems  with  flows  of  the  form 

x(k  + 1)  =  Aix{k)  +  Bi  +  UiS{k ),  (4) 

with  Ui  E  7 Znxq  and  S  E  loo-  Both  the  necessary  condition  and  the  sufficient  condition  derived 
in  this  section  hold  for  this  type  of  uncertain  system.  Robust  simulation  for  noisy  systems 
calculates  all  possible  trajectories  for  all  possible  noise  signals.  The  LMI  based  stability 
conditions  are  not  applicable  for  this  class  of  systems. 

1.4.3  Examples 

Two  simple  2-dimensional  examples  show  the  limitations  of  the  LMI  approach  and  the  value 
of  the  robust  simulation  based  algorithms.  Consider  two  piecewise  linear  systems,  each  with 
two  regions  (1  =  2)  whose  flows  are  given  by 

=  f  0.81  0.22  1  [  0.81  “0.22  ' 

“  [  -0.22  1.17  J  ’  Al  ~  [  0.22  1.17  _ 
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Figure  2:  Stable  PL  system 


with  Bq  =  B\  =  [0;0].  The  piecewise  linear  systems  are  defined  on  the  rectangle  |mi  |  <  100 
and  \x2\  <  100.  For  both  systems,  the  initial  region,  <S0,  is  the  rectangle  defined  by  |xi|  <  10 
and  \x%\  <  10  and  the  final  region,  St,  is  the  rectangle  defined  by  x,x\  <  0.1  and  x2  <  0.1.  It 
is  important  that  the  initial  region,  <S0,  is  not  the  entire  region  of  definition  of  the  PL  system. 
In  both  of  these  examples,  some  trajectories  leave  <So  and  later  reenter  it.  If  the  trajectory 
exited  the  system’s  region  of  definition,  it  would  be  unstable,  by  definition.  The  number  of 
time  steps,  t,  for  both  simulations  is  100. 

Since  Bo  =  B\  =  [0;  0] ,  the  LMI  approach  is  applicable  to  any  piecewise  linear  system 
with  these  flows.  However,  both  (2)  and  (3)  have  no  feasible  solution  and  no  conclusion  can 
be  drawn.  Though  both  Ao  and  A\  are  separately  stable,  their  combination  may  not  be. 

The  first  system  follows  the  switching  law 


x[k  +  1] 


Aox[k\  X2  >  0 
Aix[fc]  X2  <  0 


A  sample  trajectory,  shown  in  figure  2,  converges  to  the  origin,  staying  close  to  the  switching 
surface.  The  system  satisfies  the  robust  simulation  based  sufficient  condition  and  is  thus 
verified  to  be  stable.  The  necessary  condition  is  also  satisfied,  as  expected. 

The  second  system  follows  the  switching  law 


x 


A0x[k] 

Aix[A;] 


x\  >  0 
x\  <  0 


A  sample  trajectory,  shown  in  figure  3,  diverges,  oscillating  along  the  switching  surface. 
The  system  fails  the  robust  simulation  based  necessary  condition,  and  is  thus  verified  to  be 
unstable.  The  systems  also  fails  the  sufficient  condition. 
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Figure  3:  Unstable  PL  system 


In  both  of  these  examples,  the  switching  law  was  critical  to  the  stability  of  the  system.  In 
any  system  where  the  switching  law  matters,  the  LMI  stability  conditions  will  never  give  any 
information  since  it  explicitly  ignores  the  switches.  Though  2-dimensional  systems  can  be 
analyzed  by  a  variety  of  other  methods,  those  methods  do  not  generalize  to  higher  dimension. 
Both  the  robust  simulation  based  algorithms  and  the  LMI  approach  generalize  to  arbitrary 
dimension. 

A  third  2-dimensional  example  demonstrates  how  bounded  noise  can  drive  a  piecewise 
linear  system  unstable.  Consider  the  piecewise  linear  system 

(  Aox[k]  |mi  |  <0.1 

x[k  +  1]  =  <  Aixf/c]  +  B\  +  U\5[k ]  aq  <  —0.1 
\  A.2x[/c]  +  #2  mi  >  0.1 


where 

1  1  ' 

1  1  ’ 


The  piecewise  linear  system  is  defined  on  the  rectangle  |mi  |  <  100  and  \x2  <  100.  For  the 
stability  question,  the  initial  region,  So,  is  the  rectangle  defined  by  |aq|  <  10  and  |m2 j  <  10. 
and  the  final  region,  St,  is  the  rectangle  defined  by  |xi|  <  0.1  and  |m2 1  <  0.1.  The  number  of 
time  steps,  t,  is  100.  Since  some  flows  have  affine  (13,)  terms,  the  LMI  based  approach  is  not 
applicable. 


Aq  = 


0.9 

0 


Bx  = 


0 

0.9 

0.1 

0.1 


Ai  =  A2  = 


,  b2  = 


-o.: 

o.i 
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With  no  noise  (U\  =  [0;0]),  the  system  satisfies  the  robust  simulation  based  sufficient 
condition  and  is  stable.  As  shown  in  figure  4,  flows  starting  with  |xi|  >0.1  move  towards 
x\  —  0,  and  then  converge  to  the  origin.  In  fact,  once  |®i|  <0.1,  the  state  converges  to  the 
origin  quadratically. 

When  bounded  noise  is  added  as  U\  =  [0.1;  0],  the  system  is  destabilized.  By  choosing  an 
initial  condition  in  region  1  and  a  noise  signal  <5[&]  =  —1,  the  system  diverges,  as  shown  in 
figure  5.  This  system  fails  both  the  sufficient  condition  and  the  necessary  condition,  and  we 
conclude  that  it  is  unstable. 

The  next  example  demonstrated  difference  between  nominal  stability  and  robust  stability 
for  PL  systems.  Though  all  initial  conditions  converge  to  the  origin  in  the  absence  of  noise,  a 
bounded  noise  can  destabilize  the  system.  No  other  techniques  exist  for  analyzing  PL  systems 
with  noise.  Robust  simulation  accounts  for  all  possible  noise  signals  and  provides  the  answer 
in  reasonable  time. 

To  demonstrate  the  algorithms  on  larger  systems,  a  5th  order  continuous  linear  system 
with  an  input  saturation  was  randomly  generated  and  discretized.  Since  the  system  has  a 
saturation,  some  regions  have  affine  terms  and  the  LPV  methods  are  not  applicable.  Rather 
than  examine  stability  directly,  the  accuracy  of  robust  simulation  was  examined. 

The  simulation  started  with  the  hypercube  |a:[0]|oo  <  1-3  and  ran  for  20  time  steps.  At 
each  time  step,  the  difference  between  robust  simulation  and  the  exact  answer  was  measured. 
The  relative  error,  defined  as  {\xSimuiated\oo  ~  Inexact  loo) /Inexact  I  oo)>  is  shown  in  figure  6. 

For  this  example,  the  difference  between  the  exact  solution  and  robust  simulation  is  neg¬ 
ligible.  While  there  are  systems  for  which  basic  robust  simulation  differs  greatly  from  the 
exact  answer,  most  discretizations  of  continuous  systems  give  good  results.  While  this  is  a 
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relatively  small  system,  it  gives  hope  that  the  algorithms  will  perform  well  on  large  systems. 

1.5  Nonlinear  MPC  Lower  Bound 

Model  predictive  control  (MPC)  for  nonlinear  systems  typically  involves  a  non-convex  op¬ 
timization  problem.  As  with  all  non-convex  optimizations,  a  local  minimum  is  found,  but 
nothing  can  be  said  about  the  global  minimum.  With  a  careful  choice  of  cost,  constraints, 
and  system  representation,  robust  simulation  gives  a  lower  bound  on  the  optimal  cost.  If 
this  bound  differs  greatly  from  the  MPC  cost,  then  additional  optimization  may  be  desired. 
Furthermore,  the  robust  simulation  results  can  be  used  to  initialize  additional  MPC  optimiza¬ 
tions.  This  technique  is  demonstrated  on  a  simple  example. 

1.5.1  Model  Predictive  Control 

Model  predictive  control  (MPC),  also  known  as  receding  horizon  control,  is  a  technique  where 
an  on-line,  open- loop  control  problem  is  solved  at  each  time  step. [8]  Using  the  current  state, 
an  input  sequence  is  calculated  to  minimize  a  cost  while  satisfying  specified  constraints.  Only 
the  first  element  of  the  sequence  is  used,  and  then  the  algorithm  is  applied  again,  beginning  at 
the  new  current  state.  For  general  nonlinear  systems,  this  technique  results  in  a  constrained 
non-convex  optimization.  [21]  At  each  time  step,  a  local  minimum  of  the  cost  is  found,  but  no 
information  about  the  global  minimum  is  provided. 

By  restricting  the  class  of  systems  considered,  iterative  robust  simulation  can  be  used  to 
find  a  lower  bound  on  the  global  minimum.  If  the  MPC  cost  and  the  robust  simulation  bound 
are  similar,  then  additional  optimization  will  yield  little  benefit.  If  they  differ  greatly,  then 
additional  optimization  may  yield  large  improvements.  Robust  simulation  results  are  used  to 
generate  initial  conditions  for  additional  optimization  runs,  which  will  generally  converge  to 
a  better  local  minimum. 

To  allow  for  analysis  by  robust  simulation,  we  require  that  the  system  is  piecewise  linear 
(PL).  PL  systems,  described  in  more  detail  in  [15]  and  [34],  behave  as  different  affine  systems 
in  various  regions  of  state  space.  In  any  region  i  of  state  space,  the  state  transition  map  is 
defined  by 

x[k  +  1]  =  Aix[k\  +  Ai  +  Biu[k\.  (5) 

Secondly,  we  require  the  input  constraint  \u[k]\  <  1  V&.  Finally,  the  cost,  J,  must  be  a 
weighted  infinity  norm  on  the  state  error  trajectory.  Specifically, 

J=  max  \W[k](x[k]  —  xolkjjloo.  (6) 

k  e  /co . .  .k  y 

Through  W[k\,  we  can  penalize  later  errors  more  than  the  initial  transient  response  and 
penalize  certain  directions  in  state  space  more  heavily  than  others.  Level  sets  of  this  cost 
define  a  tube  through  which  the  state  trajectory  must  pass. 
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1.5.2  Algorithm 

The  idea  of  the  algorithm  is  to  find  the  largest  cost  J,  denoted  Jopt ,  such  that  there  are  no 
feasible  trajectories  with  cost  less  than  J.  This  is  accomplished  by  iterative  robust  simulation. 
Robust  simulation,  defined  in  [15]  and  described  in  more  detail  in  [12],  allows  the  calculation 
of  all  feasible  trajectories. 

The  algorithm  uses  the  following  notation: 

0  ^  Jib  —  Jopt  5:  Jmpc 

Jmpc  is  the  cost  from  the  MPC  optimization  and  represents  the  lowest  cost  from  a  single, 
known  feasible  trajectory.  Jopt  is  the  unknown  optimal  cost.  J«,  is  a  lower  bound  on  Jopt . 
The  main  goal  of  the  algorithm  is  to  find  a  large  Jw  with  a  secondary  goal  being  to  reduce 
the  gap  between  and  Jmpc. 

The  first  step  is  to  construct  an  MPC  control  sequence  and  calculate  its  cost,  Jmpc.  Since 
this  is  a  feasible  trajectory,  it  provides  an  upper  bound  on  Jopt.  Initially,  Jib  =  0,  since  0  is 
always  a  lower  bound. 

The  second  step  is  to  choose  a  cost  Jraom,  Jib  <  Jnom  <  Jmpc  and  determine  if  any 
feasible  trajectories  meeting  that  cost  exist.  This  is  accomplished  using  a  modified  robust 
simulation  algorithm  that  calculates  the  set  of  states  that  can  be  reached  at  each  time  step 
while  satisfying  the  cost  constraint. 

Robust  simulation  is  a  discrete  simulation  technique  that  calculates  all  reachable  states 
at  time  ki+\  from  a  set  of  states  defined  at  time  kt.  When  finding  MPC  lower  bounds,  we  are 
only  concerned  with  states  that  are  reachable  and  meet  the  cost  constraint.  With  this  added 
requirement,  robust  simulation  is  applied  to  the  intersection  of  the  set  of  states  reachable 
at  time  ki  and  the  set  of  states  meeting  the  cost  constraint.  Forming  this  intersection  is 
equivalent  to  finding  the  intersection  of  two  polyhedra  and  can  be  written  as  a  linear  program. 
When  this  intersection  is  empty,  there  are  no  trajectories  satisfying  J  <  Jnom  and  the  robust 
simulation  algorithm  is  terminated. 

If  no  feasible  trajectory  satisfying  J  <  Jnom  exists,  then  Jnom  is  a  lower  bound  on  Jopt. 
If  Jnom  is  an  acceptable  lower  bound  (Jnom  «  Jmpc )  then  we  are  done.  If  Jn0m  is  too  small, 
then  set  Jib  —  Jnom  and  repeat  the  second  step. 

If  robust  simulation  generates  potential  trajectories  that  meet  the  cost  constraint,  then 
two  options  exist.  Either  the  results  from  the  robust  simulation  can  be  used  to  find  a  better 
MPC  sequence  (reduce  Jmpc)  or  the  robust  simulation  step  can  be  repeated  with  a  smaller 
Jnom • 

The  algorithm  terminates  when  one  is  satisfied  with  J^.  Typically,  this  is  when  Jmpc  — 
is  small.  Several  decisions  in  the  algorithm,  including  the  choice  of  Jnomi  are  heuristic  based. 
The  use  of  robust  simulation  results  to  initialize  added  optimizations  involves  using  branch 
and  bound  on  the  input  space. 

Since  robust  simulation’s  computational  cost  grows  polynomially  with  problem  size,  this 
algorithm  also  exhibits  polynomial  growth.  Other  techniques  for  finding  better  local  minima, 
such  as  gridding  the  parameter  space,  have  exponential  computational  growth. 
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Figure  7:  cost  J  as  a  function  of  r/[l]  and  u[2] 


1.5.3  Example 


Consider  the  piecewise  linear  system 


f  rr [A;]  +  0.15u[£;],  x[k]  <  1 
\  x[k]  +  u[k],  x[k]  >  1 


with  input  constraints  |u[A]|  <  1.  Starting  from  a;[0]  =  0.95,  we  want  to  minimize  the  cost 
J  =  max{0.1a;[l],  m [2] } .  The  cost,  J,  as  a  function  of  the  control  inputs  u[l]  and  u[2]  is  shown 
in  Figure  7. 

Any  optimization  scheme  using  a  gradient  descent  algorithm  and  starting  from  u[l]  <1/3 
will  fall  into  the  local  minimum  at  u[l]  =  u[ 2]  =  -1.  MPC  based  control  schemes  are 
vulnerable  to  failures  of  this  sort  in  their  optimization  step  [24] . 

The  local  minimum  at  u[l]  =  it[2]  =  -1  has  a  cost  J  =  0.65.  This  is  more  than  a  factor  of 
6  greater  than  the  global  minimum  cost  of  0.1,  which  occurs  at  w[l]  =  1/3,  —1  <  u[2\  <  —0.9. 

Starting  the  MPC  optimization  from  u[l]  =  u[2]  =  0  converges  to  the  local  minimum, 
resulting  in  Jmpc  —  0.65.  Iterative  robust  simulation  gives  a  lower  bound  of  =  0.099.  The 
gap  between  Jmpc  and  J/b  indicates  that  the  MPC  solution  may  only  be  a  local  minimum.  A 
single  branch  and  bound  step  of  the  robust  simulation  suggests  that  u[l]  =  0.5,  ^[2]  =  0  is  a 
better  initial  condition  for  the  optimization.  A  gradient  descent  algorithm  applied  from  this 
point  quickly  converges  to  a  global  minimum,  giving  a  new  cost  Jmpc  =  0.1.  Since  «  Jmpc, 
further  optimizations  will  yield  little  benefit. 

While  this  is  a  simple  optimization  problem  that  is  easily  solved  by  gridding  the  parameter 
space,  it  demonstrates  the  feasibility  of  the  algorithm.  Future  work  will  include  applying  this 
method  to  more  realistic  problems. 
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2  Hierarchical  Modelling 

For  modeling  complex  systems,  it  is  natural  to  reduce  the  system  into  subsystems  and  model 
each  subsystem.  The  approach  taken  in  this  work  is  that  it  is  desired  for  a  model  to  be 
consistent  with  the  modeling  methodology.  Further  it  is  important  to  explicitly  represent  the 
inaccuracies  of  the  model  as  part  of  the  model. 

A  hierarchy,  interconnection  structure,  and  a  fundamental  component  data  type  are  pro¬ 
posed  and  the  choices  motivated.  The  framework  is  proposed  with  the  intention  of  being 
implemented  on  a  computer  and  having  a  family  of  models  of  different  resolution  represent¬ 
ing  a  system. 

There  are  many  issues  with  this  problem.  One  of  the  issues  is  variable  granularity  of 
the  hierarchical  model.  What  is  meant  by  granularity  is  inherently  scale.  Do  I  have  a  simple 
aggregate  model  for  a  complex  system,  or  do  I  need  to  model  each  of  it’s  components  in  detail? 
With  this  sort  of  variable  resolution,  there  are  associated  changes  in  the  state  dimension  of  the 
model  of  the  system.  These  sorts  of  tradeoff  in  choosing  a  model  should  occur  in  a  distributed 
manner. 

A  more  specific  and  well  defined  issue,  is  the  hierarchical  modelling  of  static  nonlinearities 
in  a  dynamical  system.  With  a  notion  of  the  complexity  for  the  approximation  of  a  nonlinear¬ 
ity,  there  is  a  tradeoff  between  fidelity  and  complexity.  One  implementation  of  this  solution 
is  to  model  using  sets  of  uncertain  piecewise  linear  (PL)  systems  and  simulate  using  robust 
simulation.  Sets  of  uncertain  PL  models  are  simple  to  create  and  admit  a  simple  partial 
ordering  of  their  accuracy.  Robust  simulation  explicitly  accounts  for  model  uncertainty  and 
allows  for  switching  between  models  during  a  single  simulation.  These  techniques  present  a 
practical,  computationally  tractable  solution  to  the  problems  encountered  when  simulating 
large  systems. 

2.1  Background 

2.1.1  Linear  Fractional  Transformations 

A  linear  fractional  transformation  is  shown  in  Figure  8  for  the  map  y  =  (A  *  M)u  where 
A —  D  +  CA(I  —  AA)”1!?.  In  general,  A  represents  uncertainty  and  dynamic  elements, 


Figure  8:  Linear  Fractional  Transformation 


and  M  = 


A  B 
C  D 


is  a  realization  of  the  map  A  *  M.  The  LFT  framework  results  in  a 


convenient  method  for  adding  various  types  of  uncertainty  to  systems  [42]. 
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If  A  =  f  operator  then  the  system  shown  in  Figure  8  is  just  the  standard  state  space 
representation.  LFT  systems  are  a  natural  generalization  of  state  space  representations.  By 
allowing  A  to  represent  more  general  systems  operators,  LFT  systems  provide  a  convenient 
framework  to  add  various  types  of  uncertainty  operators,  like  Nonlinear,  LTV,  LTI,  LPV,  etc 
[28]  [40]  [6],  in  which  essentially  all  the  major  state  space  results  can  be  generalized  [1]. 


2.1.2  Implicit  LFT  Systems 

An  implicit  LFT  system  is  described  by  0  =  (A *M)w  as  shown  in  Figure  9,  where  w  contains 
all  the  system  variables  ie.  there  is  no  distinction  between  inputs  and  outputs.  Implicit  LFT 
systems  are  a  generalization  of  the  behavioral  framework  proposed  by  Willems  [37]. 


Figure  9:  Implicit  LFT  system 


Any  LFT  system  can  be  converted  into  an  implicit  LFT  system.  In  first  principles  mod¬ 
eling,  like  F  =  Ma,  neither  F  or  a  is  assumed  to  be  an  input  or  an  output.  Once  either  F  or 
a  is  defined,  both  are  defined.  This  is  natural  within  the  implicit  LFT  form  but  doesn’t  fit 
and  makes  interconnections  difficult  within  the  LFT  form. 

For  tree  structured  hierarchical  models,  at  the  interconnections  there  isn’t  a  notion  of 
signal  flow.  The  interconnection  variables  become  internal  variables  to  the  system  rather 
than  inputs  or  outputs.  So  it  is  more  natural  to  not  make  the  distinction  between  inputs  and 
outputs  in  modeling. 

Within  the  implicit  LFT  framework,  the  interconnection  of  systems  is  simple.  The  implicit 
description  of  a  system  describes  the  equations  a  system  must  satisfy  (0  =  (A where 


Mi 


Ai  Bi 

Ci  D, 
and  in  addition  t 


).  So  if  two  systems  are  connected  then  they  still  satisfy  the  same  equations, 

tie  interconnection  must  be  defined  (T\W\  +  T2W2  =  0)  which  defines  the 
intersection  of  behaviors  [37].  So  the  implicit  LFT  model  of  the  interconnected  system  is 
given  by  0  =  (A  *  M)w  where: 


‘  Ai 

0 

Bi 

0  ‘ 

0 

A2 

0 

b2 

M  = 

Ci 

0 

Di 

0 

0 

c2 

0 

d2 

0 

0 

Ti 

t2 

(7) 


A  = 


Ai 

0 


0 

A2 


(8) 
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Implicit  representations  can  be  interconnected,  manipulated  and  reduced  without  committing 
to  a  particular  input-output  form,  which  is  only  relevant  to  certain  applications,  and  which 
can  be  easily  derived  a  posteriori  if  necessary  as  shown  in  D ’Andrea  and  Paganini  [5]. 

The  implicit  LFT  framework  also  allows  our  model  to  include  integral  quadratic  con¬ 
straints  (IQCs)[22j.  IQCs  are  inequalities  involving  a  quadratic  form  in  signal  space: 

/OO 

w(uj)*Hw(u)duj  <  0  (9) 

-oo 

where  II*  is  an  LTI  operator.  IQCs  can  be  used  to  define  sets  in  which  signals  must  exist, 
like  defining  sets  of  allowable  noise.  Given  n(eja;)  =  €  L ^  ^  3k  >  0  3  kl  +  II  >  0. 

By  doing  a  spectral  factorization  [42]  of  kl  +  II,  we  find  a  Q  3  II  =  kl  —  Q*Q.  Defining 
P  =  kl!2I  ^  II  =  P*P  -  Q*Q.  Equation  9  becomes  ||  Qw  ||2>||  Pw  ||2-  This  can  be  written 
as  an  uncertain  implicit  equation  of  the  form  (P  +  A cQ)w  =  0  where  Ac  is  an  arbitrary 
contractive  operator  [29]. 

2.2  General  Paradigm 

A  modeling  framework  involves  choices.  This  work  presents  and  motivates  the  choices  made 
in  this  uncertain  hierarchical  modeling  framework.  The  motivation  for  choosing  hierarchical 
modeling  is  that  not  all  components  of  a  system  are  equally  significant  and  this  should  be 
reflected  in  the  model.  In  the  case  of  a  car,  the  ashtray  is  a  less  significant  component  than 
the  engine.  If  a  hierarchical  model  is  not  used  then  all  the  system  equations  are  written  at 
the  same  level,  and  it  may  be  difficult  to  immediately  distinguish  the  dominant  dynamics  of 
the  system  from  the  trivial  dynamics.  The  hierarchy  hopefully  provides  a  quick  and  efficient 
method  for  identifying  the  significance  of  dynamics  and  doing  model  reduction  as  desired. 

The  second  choice  made  is  doing  uncertain  modeling.  For  a  real  system,  finding  an  exact 
model  is  impossible  and  characterizing  the  inexactness  is  critical  for  making  guarantees  on 
system  performance.  In  the  case  of  a  resistor,  there  is  uncertainty  on  the  exact  value  of 
the  resistance  and  the  parasitics.  When  a  model  is  reduced,  the  reduced  dynamics  must 
be  covered  with  uncertainty.  This  may  be  useful  when  cruder  models  with  uncertainty  are 
sufficient  for  a  particular  application  which  may  lead  to  reduced  computation. 

The  next  choice  made  is  a  tree  structured  hierarchy.  The  motivation  for  a  tree  structure  is 
that  a  complicated  system  is  naturally  decomposed  into  an  interconnection  of  simpler,  more 
tractable  subsystems  and  each  subsystem  can  be  similarly  reduced.  This  self  similarity  leads 
to  a  tree  structure.  For  example,  modeling  all  the  facets  of  a  car  is  quite  a  task,  but  a  more 
natural  approach  is  to  break  up  the  car  into  more  tractable  components  and  model  them 
individually.  So  rather  than  modeling  the  entire  car,  a  car  is  an  interconnection  of  an  engine, 
transmission,  exhaust  system,  cooling  system,  suspension,  etc  and  each  of  the  components  is 
modeled  separately. 

A  benefit  of  the  tree  structure  is  the  connection  with  the  object  oriented  philosophy.  The 
tree  defines  how  the  components  interact  with  each  other.  Which  simplifies  future  modifica¬ 
tions,  because  if  a  system  is  modified  only  the  modified  components  needs  to  be  remodeled, 
and  the  other  component  models  will  remain  intact.  So  in  the  car  example,  if  the  engine  is 
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replaced  with  a  different  model,  the  entire  system  doesn’t  have  to  be  remodeled.  The  old 
engine  model  is  replaced  by  the  model  of  the  new  engine. 

Another  choice  that  is  made  is  that  the  model  will  be  constructed,  implemented  and 
used  on  a  computer.  As  a  result,  the  data  structures  should  be  convenient  and  tractable 
for  computer  implementation  as  opposed  to  writing  them  out  by  hand.  The  reason  for  this 
choice  is  that  as  more  complicated  systems  have  more  and  more  detailed  models  and  are  to 
be  analyzed  or  simulated  it  is  intractable  to  do  them  any  other  way. 

The  proposed  framework  is  defined  by  a  hierarchical  tree  structure  of  the  components,  an 
interconnection  structure  of  the  components,  and  a  fundamental  data  type  for  a  component. 
Throughout  this  work,  an  inductor  is  used  to  demonstrate  the  features  of  the  framework.  An 
inductor  is  a  simple  example,  but  is  necessary  to  make  the  presentation  tractable. 

2.2.1  Component  Modeling 

Consider  the  ideal  inductor  shown  in  Figure  10  with  the  equation  d/dt(Li)  =  v  with  in- 

I 

+ 

V 


Figure  10:  Inductor 


ductance  L ,  where  v  is  an  input, 
resulting  LFT  model 


is  an  output,  and  ([>  is  the  flux  in  the  inductor. 


'  A 

B  ' 

0 

1 

C 

D 

1/L 

0 

The 

and 


Figure  8)  of  the  system  is  given  by  u  =  v,  x  =  (f),  y  =  i,  A  =  f 

As  a  result  of  choosing  this  model  of  an  inductor,  two  problems 


arise.  First,  if  the  effects  of  uncertainty  in  L,  nonlinearities  and  parasitics,  are  to  be  investi¬ 
gated,  there  is  no  convenient  way  to  do  this  without  starting  over.  Second,  by  assuming  that 
v  is  an  input  and  i  is  an  output,  our  model  may  be  incompatible  with  other  components  with 
which  it  is  interconnected  as  shown  in  Figure  11.  Figure  11  represents  an  input /output  inter¬ 
pretation  of  two  inductors  connected  in  series  («i  =  i2).  The  connection  of  the  two  outputs 
is  not  properly  defined  within  the  input/output  framework. 

The  solution  to  the  first  problem  is  easily  addressed  in  the  LFT  framework.  For  the 
inductor  example  in  Figure  10,  the  modeled  equations  should  be  replaced  by  d(j)/dt  =  v  and 
4>  —  L(I)  =  Lo(l  +  S)i.  Where  Lq  is  a  constant  which  represents  the  nominal  value  of  the 
inductor  and  S  is  an  unknown  operator.  The  model  of  the  uncertain  inductor  is  given  by 
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Figure  11:  Series  Interconnection  of  Inductors 


y  =  *,  X  =  [<f>,X2\,  U  =  V, 

'  A 

B  ' 
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T— 1 
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-1/To 
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. c 

D 
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1/Lo 
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not  a  state  in  the  conventional  sense  but  the  output  5  [1].  5  could  represent  time  variation 
in  the  inductance  due  to  saturation,  variations  in  the  core  geometry,  magnetic  links  to  other 
components,  nonlinearities,  etc. 


LFTs  provide  a  flexible  modeling  framework  for  incorporating  uncertainty  descriptions, 
but  the  input-output  assumption  is  not  desirable  for  the  modeling  of  interconnected  systems. 
It  is  not  a  priori  known  which  variables  should  be  treated  as  inputs  and  which  should  be 
treated  as  outputs  [37]. 

To  address  the  second  problem,  systems  will  be  represented  in  an  implicit  LFT  form  where 
0  =  (A'kM)w  as  shown  in  Figure  9.  There  is  no  issue  of  compatibility  between  implicit  LFT 
models  of  components  because  no  input-output  partition  has  been  made. 

For  the  inductor  example,  the  implicit  model  of  the  inductor  is  given  by  d(j)/dt  —  u, 
Lo{  1  +  £)(*)  =(f>,w  =  [«,*], 
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,  and  A  = 


A 

dt 
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o 
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For  the  fundamental  components  in  this  modeling  system,  the  model  must  be  general 
enough  to  describe  any  situation  which  may  occur.  The  model  should  have  some  information 
about  the  component,  so  that  each  component  isn’t  just  an  arbitrary  operator.  The  informa¬ 
tion  is  contained  in  the  nominal  value  of  the  component.  No  irreversible  assumptions  should 


20 


Doyle 


F49620-94-1-0420 


be  made.  For  example,  the  standard  circuit  equation  for  an  inductor  is  V  =  L  •  di/dt.  If 
this  were  the  fundamental  model  of  an  inductor,  the  assumption  has  been  made  that  L  is  a 
constant.  If  L  were  blindly  replaced  by  L(t)  the  resulting  model  would  be  incorrect.  The 
actual  equations  for  an  inductor  are  given  in  equations  10  and  11. 


*<s> 

II 

(10) 

-to 

II 

(11) 

Once  we  have  our  model  and  are  ready  to  do  analysis,  synthesis,  or  simulation  any  assumptions 
can  be  made  about  our  uncertainty  like  Al  is  a  real  parameter,  a  bounded  operator,  etc,  but 
this  is  after  the  modeling  process  and  a  part  of  the  identification  process. 

We  want  to  make  the  weakest  possible  assumptions  about  our  inductor  (L,&),  so  that  wide 
variety  of  uncertainty  assumptions  can  be  made  at  the  analysis  level.  L  and  k  are  assumed 
to  be  a  non-commuting  indeterminates  (NCIs),  ie.  it  could  be  an  arbitrary  nonlinear,  time- 
varying  operator.  For  a  real  inductor,  it’s  nominal  value  is  of  use  in  describing  it’s  operation. 
For  modeling,  L  :=  Lo  +  Al  and  k  1  +  Av,  where  A l  and  Av  are  NCIs  and  Lo  is  a  “place 
holder”  for  the  nominal  inductance.  Lo  is  used  to  describe  the  ideal  model  of  an  inductor. 
NCI’s  act  as  “place  holders”  for  uncertainty  descriptions.  It  is  important  to  note  that  by 
setting  L  Lo  +  Al  we  have  not  committed  to  anything.  This  can  be  undone  by  defining 
Al  :=  — Lq  +  A Lnew  The  Lo  term  is  added  because  the  nominal  value  is  presumably  useful 
in  describing  the  operation  of  the  system. 

Part  of  the  modeling  process  is  the  addition  on  new  components  (like  adding  parasitics 
to  the  model  of  a  inductor).  It  is  desired  that  the  current  model  can  be  refined  to  arrive 
at  a  new  model  rather  than  discarding  the  model  and  starting  over.  As  a  result  the  model 
must  have  “place  holders”  for  defining  new  interconnections.  If  these  “place  holders”  do  not 
exist  then  it  is  possible  that  at  an  interconnection  an  incorrect  equation  may  result  as  shown 
in  Figure  12.  Figure  5a  would  lead  to  the  constraint  that  i\  —  but  if  later  an  additional 
component  is  added  as  in  Figure  5b  then  i\  =  i2  +  These  two  constraints  lead  to  23  —  0 
which  may  be  incorrect. 

In  the  circuit  case,  these  “place  holders”  would  correspond  to  error  currents  for  each  node 
and  error  voltages  for  each  branch  connecting  nodes.  These  are  needed  to  describe  all  the 
possible  interconnection  of  circuits.  When  we  do  the  analysis  these  error  place  holders  must 
be  resolved.  They  can  be  set  to  0,  ie.  no  additional  current  into  a  node,  they  can  be  considered 
external  noise  terms  and  constrained  according  to  IQCs,  or  used  to  cover  unmodeled  dynamics. 

2.2.2  Component  Data  Type 

A  component  is  chosen  to  be  modeled  by  an  implicit  LFT  system  (Figure  9)  with  special 
structure.  The  implicit  LFT  structure  is  chosen  because  it  provides  a  powerful  framework  for 
describing  uncertainty  plus  it  is  natural  for  interconnecting  systems.  The  system  variables  w 
are  partitioned  as: 

^  =  [  Wn  lm  ]  (12) 
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Figure  12:  Incorrect  and  Irreversible  interconnections 


Where  wm  are  the  manifest  variables  of  the  component.  They  represent  the  variables  that  are 
connected  to  the  parent  component.  For  an  inductor,  wm  =  [  v  i  ].  wn  are  the  nominal 
manifest  variables,  and  lm  are  the  latent  manifest  variables.  wn  are  used  to  describe  the 
nominal  dynamics  of  the  component.  wn  and  lm  are  used  to  describe  the  interconnection  of 
the  nominal  component  with  the  child  components.  lc  are  the  latent  crosstalk  variables,  and 
are  the  terms  used  to  describe  crosstalk  connections  (connection  between  unrelated  compo¬ 
nents).  le  are  the  latent  error  variables,  they  are  used  to  describe  interconnection  errors  and 
act  as  “place  holders”  for  describing  future  interconnections.  It  is  difficult  to  motivate  this 
partition  of  the  system  variables,  but  this  partition  should  be  clear  after  the  hierarchy  and 
interconnection  structure  of  the  system  is  defined. 

The  special  structure  of  A  *  M  is  described  by: 

B  =  [  0  h  0  b2  63  ]  (13) 


0  d\  0  g?2  0 
dm  dn  d[  dc  de 


A  =  diag 


’  7  Qnl)  Ai,  •  •  • ,  Ar 


(14) 

(15) 

(16) 


A,£?,C,and  D  are  constant  matrices,  qi  are  arbitrary  constants.  A;  are  noncommuting  inde- 
terminates.  A  is  used  to  describe  the  dynamics  and  uncertainty  of  the  nominal  component. 
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Within  A  the  standard  dynamic  element  is  Jr  instead  of  the  standard  f.  The  reason  for  this 
choice  is  that  f  is  not  an  operator  because  an  initial  condition  must  be  specified. 

The  component  model  of  a  “real”  inductor  is  shown  in  Figure  13.  The  implicit  LFT 


Figure  13:  A  Model  of  a  “Real”  Inductor 


model  is  given  by  (using  Figure  9, equations  13,  14,  15, and  16:  wm  =  [v,i},wn  =  [vl,  *z]> 
lm  =  [vR,iR,vi],  lc  =  [],  le  =  [vei,ve2,iei,ie2,<t>e\  (vei  corresponds  to  a  voltage  error  of  R,  ve2 
corresponds  to  a  voltage  error  of  L,  and  </>e  corresponds  to  the  flux  error  in  the  inductor), 
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The  terms  that  are  not  listed  are  empty  ie.  =  []. 
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2.3  Hierarchy  of  Components 

The  hierarchical  structure  of  a  system  modeled  using  this  framework  is  a  tree  (Figure  14).  At 
each  node  of  the  tree  there  is  a  component  of  the  form  presented  in  section  2.2.2.  Components 
can  be  ranked  hierarchically  if  they  are  on  the  same  branch  of  the  tree.  The  component 
closer  to  the  root  is  considered  hierarchically  above  the  component  closer  the  leaf  of  a  branch. 
Components  not  on  a  common  branch  are  not  related  hierarchically.  A  particular  branch 
may  be  more  significant  than  another.  For  the  inductor  model  in  Figure  13,  the  hierarchy  is 


Figure  14:  Hierarchical  Structure 


The  link  that  connects  two  nodes  on  the  tree  contains  the  information  of  how  the  la¬ 


tent  variables  of  the  higher  component  are  connected  to  the  manifest  variables  of  the  lower 
component.  This  is  a  hierarchical  interconnection.  The  link  defines  T\  and  T2  of  (7)  for  the 
interconnection.  Assuming  that  M\  is  the  lower  component  and  M2  is  the  higher  compo¬ 
nent,  T\  =  [  I  0  0  0  0  ]  and  T2  =  [  0  0  t<z  0  0  ] .  For  the  inductor  model, 
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2.3.1  Nonhierarchical  Interconnections 

Hierarchical  interconnections  are  not  the  only  interconnections.  Any  two  nodes  may  be  in¬ 
terconnected.  The  nonhierarchical  interconnections  are  called  crosstalk.  Crosstalk  is  used  to 
model  secondary  interactions  between  components  like  magnetic  interconnection  of  disjoint 
inductors,  gravitational  interaction  of  masses,  two  disjoint  flexible  structures  linked  by  air, 
etc.  The  structure  of  the  model  of  the  system  is  a  web  (Figure  16)  rather  than  a  tree,  but  the 
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Figure  16:  Web  Structure  of  Model 


web  has  an  underlying  tree  structure  from  the  hierarchy  as  shown  in  Figure  14.  The  solid  lines 
in  Figure  16  denote  hierarchical  connections  and  the  dashed  ones  the  crosstalk  connections. 

A  non- hierarchical  interconnection  link  is  actually  a  component  describing  the  dynamics  of 
the  connection  (dynamics  of  air  linking  two  disjoint  flexible  structures).  The  link  defines  how 
the  latent  crosstalk  variables  of  the  two  hierarchical  components  are  connected  to  the  crosstalk 
component.  The  link  defines  an  interconnection  of  3  components.  So  generalizing  (7)  we  have 
Ti,  T2,  and  T3  where  T*  =  [  0  0  0  U  0  ]  for  i=l  or  2  and  T3  =  [  £3  0  0  0  0  ] . 

2.3.2  Component  Refinement 

The  process  of  component  refinement  is  an  improvement  of  the  model  of  a  component.  A  more 
accurate  nominal  model  of  the  system  (without  uncertainty,  noise,  etc)  can  be  parametrized. 
As  a  result,  the  uncertainty  description  of  the  system  doesn’t  have  to  be  as  conservative. 

Component  refinement  involves  interconnecting  a  component  to  more  “parasitic”  subcom¬ 
ponents  or  describing  crosstalk  interaction  between  components,  which  models  previously  un¬ 
described  phenomenon.  The  proposed  hierarchical  modeling  framework  was  designed  with 
this  option  in  mind.  The  idea  being  that  a  model  can  describe  new  phenomenon  without 
having  to  start  over.  This  is  the  purpose  of  the  latent  error  variables  l€. 

As  a  model  is  refined  new  latent  manifest  variables  are  created  for  components.  In  the 
circuit  case  these  are  the  new  interconnection  variables,  the  voltages  associated  with  new 
nodes,  the  error  terms  for  new  nodes  and  branches,  and  refinements  of  old  error  terms. 

In  the  process  of  refining  the  components  of  a  model  (arriving  at  a  more  detailed  model) 
there  is  an  implicit  reduction  of  the  uncertainty  necessary  to  describe  the  system  dynamics 
because  the  uncertainty  is  still  represented  by  an  NCI.  As  more  phenomenon  are  modeled, 
the  inaccuracy  of  the  model  is  reduced. 

2.3.3  Model  Reduction 

Assuming  that  we  have  a  detailed  model  of  a  system  with  all  the  nominal  values  determined, 
NCIs  characterized,  and  error  variables  resolved,  it  may  be  desirable  to  view  the  system  at  dif- 
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ferent  levels  of  resolution.  There  is  a  tradeoff  between  the  complexity  and  fidelity  of  a  model. 
When  the  model  is  used  for  analysis,  synthesis  of  a  controller,  or  simulation  less  detailed 
models  may  be  sufficient  for  the  task  at  hand  which  would  reduce  the  required  computation. 
There  is  no  reason  to  use  highly  detailed  models  when  crude  models  are  sufficient. 

In  this  framework,  model  reduction  can  be  done  efficiently.  The  model  reduction  process 
involves  paring  the  hierarchical  tree.  When  a  branch  is  cut  the  dynamics  of  the  removed 
components  must  be  covered  by  uncertainty.  The  removed  branch  would  be  replaced  by 
IQCs. 


2.4  Hierarchical  Modelling  of  Static  Nonlinearities 

When  creating  a  mathematical  model  of  a  physical  system,  any  degree  of  complexity  is  possi¬ 
ble.  A  resistor  model  can  simply  be  e  =  ir  or  can  include  a  variety  of  parasitic  effects.  When 
developing  models  for  interconnection  with  other  systems,  the  degree  of  accuracy  needed  is 
not  known  a  priori.  Effects  critical  to  one  system  may  be  superfluous  to  another.  Thus,  a 
set  of  models  with  varying  degrees  of  accuracy  are  needed.  Finally,  the  set  of  models  must 
be  useful  for  computation. 

Model  hierarchies  are  an  area  of  active  research  in  the  computer  graphics  community. 
Wavelets  are  used  to  create  sets  of  models  of  3  dimensional  objects  for  on-screen  rendering. 
As  an  object  moves  toward  the  foreground  of  an  image,  more  complex  models  are  used  to 
display  finer  detail.  [32] 

Modeling  for  control  has  its  own  set  of  requirements.  Primarily,  simulation  and  analysis 
computations  must  grow  reasonable  as  state  dimension  increases.  In  computer  graphics,  all 
problems  are  in  2  or  3  dimensions.  Secondly,  some  measure  of  model  accuracy  is  needed.  Is 
a  nonlinear  2  state  model  more  accurate  than  a  linear  5  state  model?  Piecewise  linear  (PL) 
systems  satisfy  both  of  these  requirements. 

A  PL  system,  described  in  more  detail  in  [15]  and  [34],  is  a  nonlinear  system  of  the  form 

x[k  +  1]  =  +  Ai  +  B{n[k]^  x  £  i  £  1 . . .  /.  (1 7) 

where  Ri  is  one  of  a  finite  number  of  regions  in  state  space  and  \n[k]  \  <  1.  This  type  of  system 
is  a  conceptually  simple  extension  of  a  linear  system;  each  region  of  state  space  behaves  as  an 
affine  system.  PL  systems  are  also  easily  simulated  and  implemented  on  digital  computers, 
and  can  approximate  a  nonlinear  system  to  any  degree  of  accuracy. 

A  partial  ordering  of  model  accuracy  is  determined  by  the  number  of  piecewise  linear 
regions,  £,  and  the  amount  of  noise  in  the  model.  Larger  l  and  less  noise  lead  to  more 
accurate  models.  There  are  several  methods  for  measuring  the  amount  of  noise  in  the  system. 
The  simplest  is  the  max^i.../  \Bi\.  More  complex  measures  can  weight  the  Bi  or  use  other 
norms. 

To  demonstrate  the  ideas  from  the  previous  section,  a  model  library  is  created  for  a  simple 
static  nonlinearity.  The  saturation 


2 

y  =  —  arctan  it, 

7T 


(18) 
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Figure  17:  The  saturation  nonlinearity 


shown  in  figure  17,  can  be  described  by  a  variety  of  PL  systems.  The  simplest  model  is  to 
treat  it  as  one  region,  1  =  1,  with  noise, 

y  =  Ou  +  n,  |n|  <  1  (19) 


While  (19)  does  not  describe  the  behavior  of  the  nonlinearity,  it  does  describe  its  bounds. 
This  model  is  useful  when  the  saturation  has  little  effect  on  the  overall  system. 

A  more  common  approximation  for  saturations  is  a  3-segment  representation.  When 
approximating  (18)  by  a  PL  system  with  l  =  3,  some  assumptions  must  be  made.  First,  this 
approximation  will  hold  for  all  values  u.  Second,  the  model  will  be  symmetric  about  u  =  0. 
Finally,  the  maximum  size  of  the  noise  signal  will  be  minimized.  The  model 


V  = 


< 


—0.87  +  0.13n 
0.38^  +  0.13n 
0.87  +  0.13n 


u  <  -2.32 
-2.32  <  u  <  2.32 
u  >  2.32 


(20) 


shown  in  figure  18,  satisfies  these  requirements.  Note  that  the  saturated  regions,  |u|  >  2.32 
are  not  nominally  equal  to  the  saturated  values.  For  this  model,  they  are  offset  to  reduce  the 
amount  of  noise  needed  to  cover  the  true  saturation  (18).  As  shown  by  the  dotted  lines,  there 
is  always  a  noise  that  makes  the  approximate  system  equal  the  true  system. 

A  more  accurate  model  can  be  created  by  using  a  5-segment  approximation.  Using  the 
same  assumptions  as  before,  the  model,  shown  in  figure  19,  is 


y=  s 


—0.95  +  0.05n 
0.06u  —  0.55  +  0.05n 
0.49u  +  0.05n 
0.06it  +  0.55  +  0.05n 
0.95  +  0.05n 


u  <  —6.14 
-6.14  <  u  <  -1.29 
-1.29  <  u  <  1.29 
1.29  <  u  <  6.14 
u  >  6.14 


(21) 
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All  three  of  these  models  are  valid  representations  of  original  nonlinearity,  18.  The  more 
complex  models  contain  less  uncertainty  (noise),  but  may  require  more  computation  during 
simulation. 


2.4.1  Multiresolution  simulation 


Traditional  simulation  techniques  yield  one  result  for  an  initial  condition  and  noise  signal.  If 
several  models  exist  for  a  system,  and  each  yields  a  different  simulation  result,  which  one  is 
correct?  The  simulation  technique  should  give  compatible  results.  Also,  the  ability  to  change 
models  during  simulation  is  desired.  One  model  may  be  accurate  for  one  operating  region 
and  a  second  may  describe  a  different  region  in  detail. 

When  simulating  uncertain  systems,  the  result  should  be  a  set  of  all  possible  final  com 
ditions,  not  a  single  final  condition.  Since  the  system  is  uncertain,  it  is  impossible  to  know 
the  final  condition  exactly.  The  simulation  should  also  allow  a  set  of  initial  conditions,  so  the 
results  of  a  prior  simulation  can  be  used  start  future  runs.  More  accurate  (less  uncertain) 
models  should  yield  smaller  sets  of  final  conditions  than  coarser  models  for  the  same  set  of  ini¬ 
tial  conditions.  Finally,  all  models  should  be  compatible.  For  any  set  of  initial  conditions,  the 
simulation  results  from  all  models  should  have  some  points  in  common.  Robust  simulation, 
described  in  [15]  and  [12],  meets  these  requirements. 

Robust  simulation  is  the  simulation  of  sets.  For  a  given  initial  condition  set  and  uncer¬ 
tainty  description,  the  set  of  all  possible  final  conditions  are  calculated.  For  model  libraries 
that  share  the  same  state  variables  (or  have  invertible  mappings  between  different  state  de¬ 
scriptions),  any  model  can  be  used  at  any  time  step  without  recalculating  prior  results.  This 
allows  the  simulation  technique  to  choose  the  most  appropriate  model  for  a  given  initial 
condition  set,  i.e.  multiresolution  simulation. 

The  ability  to  do  dynamic  model  selection  is  a  direct  result  of  using  robust  simulation. 
The  result  at  any  time  step  can  be  used  as  the  input  for  simulating  any  model  one  time 
step  forward.  The  technique  does  not  provide  guidelines  for  model  selection;  it  only  gives 
the  ability  to  change  models.  As  long  as  all  models  correctly  describe  the  same  system,  the 
results  are  guaranteed,  by  construction  of  the  algorithm,  to  be  compatible. 

This  simulation  technique  exhibits  polynomial  computational  growth  in  all  variables.  As 
problem  size  grows,  computational  cost  grows  reasonably.  For  each  time  step  of  simulation, 
roughly  nl 2  linear  programs  must  be  solved.  This  technique  is  only  applicable  for  models  in 
the  form  of  (17)  and  is  conservative.  If  the  results  are  too  conservative,  the  simulation  can 
be  systematically  refined  until  the  exact  solution  is  obtained. 

Three  one-dimensional  systems  are  used  to  demonstrate  the  multiresolution  simulation 
ideas.  While  more  complex  examples  can  easily  be  solved,  these  simple  show  the  possible  re¬ 
sults.  Robust  simulation,  the  algorithm  used  during  simulation,  has  already  been  successfully 
tested  on  higher  order  models.  [15] 

The  first  example  demonstrates  that  in  some  circumstances,  even  the  coarsest  model  may 
give  acceptable  results.  Consider  the  stable  system 


x[k  +  1]  =  0Ax[k]  +  0.05—  arctanufA;] 

7 r 


(22) 
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time  step 

lower  bound 

exact 

upper  bound 

0 

100 

100 

100 

1 

39.95 

39.95 

40.05 

2 

15.93 

15.93 

16.07 

3 

6.32 

6.32 

6.48 

4 

2.47 

2.48 

2.64 

5 

0.94 

0.96 

1.11 

Table  1:  Simulation  results  for  coarse  model 


Saturation  Model 

min  re  [50] 

max  a;  [50] 

Model  (19) 

-280 

2200 

Model  (20) 

-0.47 

140 

Model  (21) 

-0.13 

0.13 

Exact 

6.0  e-11 

6.5  e-8 

Table  2:  Simulation  results  for  the  unstable  system 


with  the  feedback  law  u[k\  —  —  x[k\.  This  system  is  open  loop  stable,  and  the  input  has  very 
little  control  authority.  All  simulations  are  for  5  time  steps  starting  from  x[0]  =  100. 

The  first  simulation  of  (22),  whose  results  are  shown  in  table  1,  uses  (19)  to  approximate 
to  approximate  the  saturation. 

For  this  system,  even  the  simplest  model  gives  reasonable  results.  The  final  result,  0.942  < 
x[5]  <  1.107,  differs  by  only  5%  from  its  center  value  and,  as  expected,  contains  the  exact 
result  of  0.956.  Note  that  for  the  first  few  time  steps,  the  lower  bound  is,  after  rounding, 
identical  to  the  exact  value  .  This  is  because  the  worst  case  uncertainty  (noise),  which  achieves 
the  lower  bound,  is  also  the  value  needed  to  make  the  approximation  match  the  true  model. 

Running  the  same  simulation  using  the  approximation  (20)  yields  a  final  result  of  0.942  < 
x[5]  <  0.963.  This  result  is  much  tighter  because  the  saturation  only  takes  on  values  between 
1  and  0.74  in  the  saturation  region.  Finally,  using  (21)  to  approximate  the  saturation  gives 
0.953  <  x[5]  <  0.961.  As  expected,  the  model  with  the  least  uncertainty  gives  the  tightest 
bounds. 

The  second  example  demonstrates  that  radically  different  results  can  be  obtained  from 
different,  seemingly  accurate,  models.  The  nominally  unstable  system 


x[k  +  1]  —  1.1  x[k]  H — arctanuffc] 

7 r 


(23) 


is  stabilized  about  the  origin  by  the  feedback  law  u[k]  =  —x[k\.  What  values  can  ,x[50]  attain 
if  7.5  <  s[0]  <  8.5? 

Table  2  shows  results  from  the  robust  simulation  using  each  of  the  three  saturation  ap¬ 
proximations.  Since  the  the  coarsest  approximation,  model  (19),  does  not  account  for  the 
sign  of  the  saturation,  that  model  does  not  result  in  a  stable  system.  As  shown  in  table  2, 
the  simulation  values  differ  greatly  from  the  exact  values. 
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For  model  (20),  which  appears  to  approximate  the  nonlinearity  fairly  well,  some  states 
approach  the  origin  and  other  states  diverge.  Model  (21)  shows  that  all  states  in  the  range 
7.5  <  a[0]  <  8.5  converge  to  the  origin.  For  this  example,  the  3  segment  saturation  approxi¬ 
mation  is  not  sufficient;  the  5  segment  model  is  needed. 

The  results  for  models  (20)  and  (21)  indicate  a  limitation  of  the  modeling  framework.  Since 
uncertainty  can  enter  as  a  constant  offset,  the  simulation  never  converges  to  an  equilibrium 
point,  but  always  to  a  ball  around  an  equilibrium  point.  For  example,  For  model  (20),  this 
set  is  —0.47  <  x  <  0.47.  Even  if  x[0]  =  0,  the  simulation  result  would  be  |or[50]  |  <  0.47.  The 

size  of  this  set  is  a  function  of  the  dynamics  around  the  fixed  point  and  the  size  of  the  model 

uncertainty.  Model  (21),  which  has  much  smaller  uncertainty,  reaches  the  ball  | m [50]  |  <  0.13. 

A  third  example  demonstrates  multiresolution  simulation.  Consider  the  nominally  unsta¬ 
ble  system  with  a  quantized  input 

x[k  +  1]  =  1.9a:[fc]  +  g(rt[A;]) 

q(u)  =  u  —  (u  +  0.05)  mod  0.1  +  0.05.  (24) 

While  q  can  be  exactly  represented  as  a  PL  mapping,  it  has  10  PL  regions  for  every  unit 

of  state  space  modeled.  A  PL  mapping  valid  for  |rc|  <  10  requires  200  PL  regions.  A  much 
simpler  model  for  the  quantizer  is 


q(x)  =  x  +  0.05n 


(25) 


Closing  the  loop  with  unity  feedback  and  using  approximation  (25)  for  the  quantizer  gives 

x[k  +  1]  =  0.9x[A:]  +  0.05n[fc] 


Simulating  the  initial  condition  set  |a[0]|  <  100  for  1000  time  steps  gives  the  result  |a;[1000]|  < 
0.50,  the  smallest  range  attainable  for  this  amount  of  uncertainty. 

A  substantially  better  result  can  be  obtained  by  changing  the  quantizer  approximation 
midway  through  the  simulation.  Another  piecewise  linear  model  for  the  quantizer  is 


x  +  0.05n, 

a  >  .55 

0.5, 

0.45a;  <  x  <  0.55 

0.4, 

0.35a;  <  x  <  0.45 

0.3, 

0.25a;  <  x  <  0.35 

0.2, 

0.15a  <  x  <  0.25 

0.1, 

0.05a;  <  x  <  0.15 

o, 

—0.05a;  <  x  <  .05 

0.1, 

-0.15a  <  a  <  -0.05 

0.2, 

—0.25a;  <  x  <  —0.15 

0.3, 

—0.35a;  <x<  —0.25 

0.4, 

-0.45a;  <  a  <  -0.35 

0.5, 

—0.55a;  <  x  <  —0.45 

x  +  0.05n, 

a  <  -0.55 
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This  model  exactly  represents  the  quantizer  near  the  origin,  and  uses  the  same  approximation 
as  (25)  for  |z|  >  0.55. 

After  running  the  first  100  steps  of  the  simulation  using  approximation  (25),  the  range 
| m [1 00]  |  <  0.503  is  obtained,  roughly  the  smallest  region  attainable  when  using  approximation 
(25).  Since  the  region  |x[100]  <  0.503  is  valid  for  any  model  of  the  system,  it  can  be  used  to 
initialize  the  simulation  for  a  different  model.  Switching  to  approximation  (26)  for  remainder 
of  the  simulation  gives  the  result  |rc[1000]|  <  0.095.  This  range  is  less  than  one  fifth  the  size 
of  the  result  when  only  model  (25)  is  used. 

The  multiresolution  simulation  not  only  gives  a  better  result,  but  also  greatly  reduces 
computational  cost.  While  the  same  result  could  be  obtained  by  using  model  (26)  for  the 
entire  simulation,  the  computational  cost  is  much  greater.  By  using  (25)  for  the  first  100 
steps,  about  25000  fewer  linear  programs  are  solved.  Since  the  number  of  linear  programs 
solved  grows  as  0(tnl2),  this  savings  is  even  larger  for  higher  order  systems. 

3  Robustness  Analysis 

The  structured  singular  value,  /i,  introduced  by  Doyle  in  [7]  has  proven  to  be  a  useful  frame¬ 
work  for  robustness  analysis.  The  main  advantage  of  the  paradigm  is  that  many  robustness 
problems  can  be  recast  as  /x  problems  with  a  particular  uncertainty  description. 

Although  computing  /x  is  NP-hard  [4],  computationally  tractable  upper  and  lower  bounds 
exist  [39,  26,  28,  38,  41,  2].  For  some  problems,  the  gap  between  the  upper  and  lower  bounds 
may  be  large.  Strategies  for  reducing  the  gap  are  based  upon  branch  and  bound  techniques 
[27].  The  branch  and  bound  algorithms  involve  dividing  a  /x  problem  into  two  /x  problems  and 
repeating  this  process  until  the  gap  is  acceptable.  Critical  to  these  techniques  is  preventing 
exponential  growth  in  the  number  of  active  n  problems.  For  existing  branch  and  bound 
algorithms  for  /x,  the  branching  is  limited  to  axially  aligned  branches. 

An  extension  of  the  branch  and  bound  algorithm  leads  to  probabilistic  /x  [43,  19].  Proba¬ 
bilistic  fi  seems  to  be  computationally  intractable  even  for  low  dimensional  problems  (>  4), 
the  bounds  are  very  conservative  for  any  reasonable  computation  time.  Even  though  the 
branch  and  bound  and  probabilistic  /x  algorithms  are  trivial  to  parallelize,  there  axe  still 
fundamental  improvements  that  can  be  made  in  the  algorithms. 

The  focus  of  this  work  is  developing  a  richer  class  of  uncertainty  descriptions  for  which 
analogous  bounds  to  p.  can  be  computed.  Traditionally,  this  has  meant  adding  to  the  possible 
uncertainty  structures,  like  real  parameters.  This  work  takes  a  different  approach.  Given 
a  matrix  and  a  block  structure  for  the  uncertainty,  ^  defines  a  distance  to  singularity  of 
(I  —  M A)  in  the  space  of  allowable  A.  For  the  standard  /x  problem  this  distance  is  measured 
using  the  oo-norm.  It  turns  out  that  analogous  bounds  can  be  derived  if  the  distance  is 
measured  using  the  2-norm. 

With  the  development  of  the  new  set  descriptions,  it  is  possible  to  do  more  creative 
branching,  which  may  result  in  better  algorithms  for  probabilistic  /x  and  branch  and  bound 
for  /x.  There  may  also  exist  physical  problems  which  fundamentally  have  spherical  constraints 
on  the  uncertainty,  but  this  is  not  the  main  motivation  for  the  results  in  this  work. 
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3 . 1  Preliminaries 

The  notation  is  standard.  If  M  is  a  matrix,  then  MT,  M*  denote  the  transpose  and  conjugate 
transpose  matrices,  respectively.  The  n  x  n  identity  matrix  will  be  denoted  by  Jn,  I  is  used 
when  the  dimension  is  obvious.  tn  is  an  n  x  n  matrix  of  l’s.  diag[A\,  •  •  •  ,An\  is  a  block 
diagonal  matrix  with  A{  on  the  diagonal.  The  vec(M)  is  a  vector  which  is  the  columns  of  M 
stacked  into  a  vector. 

The  spectral  radius,  real  spectral  radius,  and  maximum  positive  real  eigenvalue  of  a  square 
matrix  M  are  denoted  by  p(M),  pr(M),  and  A r(M),  respectively.  The  maximum  and  min¬ 
imum  singular  values  of  a  matrix  M  are  denoted  by  a(M)  and  a(M)  respectively.  Tr(M) 
is  the  trace  of  M.  det(M)  is  the  determinant  of  M.  The  Frobenius  norm  of  a  matrix  M  is 
defined  as  \\M\\p  =  Ylij  \mij\2'  A  Hermitian  matrix  M*  =  M  G  Cixn  is  said  to  be  positive 
(semi)  definite  if  x*Mx  >  0(>  0)  for  all  nonzero  x  E 

The  Kronecker  product  of  two  matrices,  A  and  i?,  is  denoted  by  A<S>B.  The  Hadamard  (or 
Schur)  element-wise  product  of  two  matrices  A  —  [a>ij]  and  B  =  [bij]  of  the  same  dimensions 
is  defined  as  AoB  =  In  this  work  the  Hadamard  product  is  sometimes  used  to  specify 

matrix  structure.  For  example,  to  specify  that  a  matrix  is  diagonal  the  notation  I  o  D  =  D 
is  used.  Some  useful  properties  of  the  Hadamard  product  follow. 

Lemma  1  If  A  =  I  o  A  and  B  =  I  o  B  then  AXB  =  X  o  (a  *  bT )  where  A  =  diag[a\ ,  *  -  • ,  an] 
and  B  =  diag[b\ , 


Lemma  2  (Schur  Product  Theorem)  [11]  If  X%  >  (>)0  and  X  —  XioX2  then  X  >  (>)0 

The  fundamental  picture  in  the  /i-paradigm  [28]  is  shown  in  Figure  20.  Where  M  is  a 
matrix  and  A  is  a  matrix  in  a  set  of  allowable  matrices.  The  set  of  allowable  A  is  described 
by  block  diagonal  structure,  and  a  size  constraint,  a(M)  <  1.  The  fundamental  question  in 


X 


Figure  20:  Standard  Interconnected  System 


the  p-paradigm  is  are  there  any  nontrivial  solutions  to  the  loop  equations  shown  in  Figure  20. 
p(M)  answers  that  question. 
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Definition  3  The  uncertainty  structure, 

_  {ding  SrlIkrl ,  •  •  ■  ,  SrriRIkrnR  ,  Scilkcl,  ■  •  ■  ,  ScncIkcnc  >  ^/i>  ’ ' '  )  ^fnp 
Sri  €  K,  £ci  G  C,  andA fi  G  C7^  }• 

Definition  4  For  M  G  K.nxn,/UA  is  defined  as 


(27) 


JUa(M):  = 


min{a(A)  :  A  E  A,  det(/  —  M A)  =  0} 
unless  no  A  E  A  makes  I  —  M A  singular 7  in  which  case  =  0. 


(28) 


Alternatively  equation  28  can  be  written  as 

Ma(A0  =  sup  pr(MA)  =  sup  Ar(MA)  (29) 

AgA,ct(A)<1  AGA,a(A)<l 

With  the  appropriate  projection  the  uncertainty  set  in  equation  29  defines  a  hypercube  of 
dimension  hr +  nc  +  rip- 

Computing  p&(M)  is  NP-hard  but  computationally  tractable  upper  and  lower  bounds 
exist.  The  lower  bound  computation  is  based  upon  equation  29.  Choose  any  A  E  A,  then 
pr(M A)  is  a  lower  bound  to  /zA(M).  The  lower  bound  algorithms  involve  strategies  to  choose 
A  to  maximize  Xr(MA)  and  is  a  non-convex  optimization  problem.  The  upper  bound  which 
is  a  Linear  Matrix  Inequality  (LMI), 

Ma (M)  <  pA(M)  -  DJri^G{r-M*DM  +  j(GM  -  M*G)  -  72D  <  0}  (30) 

is  a  convex  optimization  problem  where, 

D  ^JDiag  [TVi?  >  ^tur i  -^cij  •>  -Gene >  ?  d fnpInfF ]  9}  nnd  (31) 

G  =  {Diag  [Gr i,  •  •  • ,  GrriR,  0,  •  ■  • ,  0] :  Gri  =  G *^}  (32) 

To  define  /z,  a  matrix  and  a  uncertainty  structure,  A,  must  be  specified.  In  this  pre¬ 
sentation,  notation  for  standard  /z  is  p^(M).  For  uncertainty  with  spherical  and  elliptical 
constraints,  ias^(M)  and  pe^(M)  are  used.  For  a  /z  problem  with  combinations  of  spherical, 
elliptical,  and  standard  constraints,  pg^(M)  is  used.  If  a  /z  problem  is  specified  without  an 
uncertainty  structure,  the  default  uncertainty  structure  is  indicated.  The  default  uncertainty 
structure  is  specified  in  definition  5. 


3.2  Spherical  /z 

For  simplicity  in  derivation  and  presentation  it  is  assumed  that  the  uncertainty  consists  of  only 
complex  scalars,  ie  no  full  blocks,  repeated  scalars,  or  real  parameters,  and  the  uncertainty 
set  is  described  by  a  Frobenius  norm  constraint  rather  than  an  induced  norm  constraint  like 
a  (A)  <  1. 
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Definition  5  The  default  uncertainty  structure  is 

As:=  {diag[5i,- ■  ■  ,Sn]:8i  G  C}.  (33) 

The  level  curves  of  ||A||/?  are  defined  by 

X>|2  =  ^,  (34) 

i~  1 

which  describes  hyperspheres  in  As  space,  hence  the  name  spherical  p  (ps). 

Definition  6  For  M  E  Rnxn;  Spherical  / 1 ,  ps  is  defined  as 

^siM)  min|||^\||F  :  A  E  As, det (I  —  MA)  =  0} 
unless  no  A  E  As  makes  I  —  MA  singular ,  in  which  case  ps(M):=  0. 

Alternatively  equation  35  can  be  written  as 

ps(M)  =  sup  pr(MA)  (36) 

A€As,||A||jr<l 

Much  like  the  standard  problem,  computing  spherical  fi  exactly  seems  to  be  intractable. 
As  a  result,  one  is  resigned  to  computable  upper  and  lower  bounds  as  in  the  standard  p  case. 

Theorem  7  fis(M )  <  ji(M)  <  y/nps(M) 

Proof.  The  statement  follows  immediately  from 

{A:  A  E  As,||A||f  <  1}  C  {A:  A  E  As,a(A)  <  1}  (37) 

{A:  A  E  As,a(A)  <  1}  C  {A:  A  E  As,  ||A||F  <  (38) 

and  equations  36  and  29. 

o 

As  a  consequence  of  Theorem  7  any  upper  bound  of  is  also  an  upper  bound  /is(M)  but 

this  would  be  a  conservative  bound  which  scales  badly  with  dimension. 


3.2.1  An  Upper  bound  to  Spherical  p 

A  /Ja(M),  which  is  not  necessarily  an  upper  bound  of  //(M),  can  be  computed  by  solving  an 
LMI  similar  to  equation  30.  The  construction  of  the  LMI  follows. 


Lemma  8  (Schur  Complement)  [3]  The  matrix  inequality , 

>0 


Q  S 
S *  R 


with  Q  =  Q*  and  R  =  R*  is  equivalent  to 

Q  -  SR-'S*  >  0 
R>  0 


(39*) 

(40) 

(41) 
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The  following  Lemma  is  used  in  a  separating  hyperplane  argument  to  construct  a  sufficient 
LMI  condition  for  there  being  no  nontrivial  solutions  to  a  quadratic  form. 


Lemma  9  If  A  >  0  and  B  >  (>)0  then  Tr(AB)  >  (>)0. 

Proof. 


Tr(AB)  =  Tr(AB^B^)  =  Tr(B^AB^) 

(42) 

o 

Bh  =  (£s)*  =4  B*ABh  >  (>)0  =►  Tr(B^AB^)  >  (>)0 

(43) 

Theorem  10 

H,{M)  <  Jis{M)  =  mf  {7:  M*{1 0  D)M  -  7 2D  <  0} 

(44) 

Proof,  x  and  y 

are  the  output  and  input  vectors  of  A,  respectively. 

i>i2<7-2 

i— 1 

(45) 

44  <5*5  <  7-2  for  6  =  [5i  62  •  •  ■  Sn]T 

(46) 

7-2  S*  1 

44  '  T  >0  Applying  Lemma  8 

<)  /  , 

(47) 

Using  Xi  =  Siyi, 

1  T  /  a.  >0  Applying  Lemma  1 

Z  Io(yy*)  J  - 

(48) 

4$  I  o  (yy*)  —  i2 xx*  >  0  Applying  Lemma  8 

(49) 

Using  y  =  Mx, 

Hs(M)  <7<^=>y  =  o;  =  0is  the  only  solution  to  equation  49 

(50) 

Using  Lemma  9 

j 

If  3D  >  0:  Tr(D{I  o  ((Mx)(Mx)*)  -  72xx*))  <  0,  Vx  G  C" 
then  equation  49  has  no  nontrivial  solutions 

(51) 

Tr(D(I  o  (yy*)  —  j2xx*))  <  0,  Vx  G  C" 

(52) 

<S4>  x*M*(I  o  D)Mx  -  7 2x*Dx  <  0,  Vx  G  C1 

(53) 

44  M*(I  o  D)M  —  j2D  <  0 

(54) 

o 
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The  main  difference  between  this  derivation  and  the  /i(M)  formulation  is  in  deriving 
the  quadratic  form  in  equation  49.  The  rest  of  the  derivation  corresponds  to  the  standard 
separating  plane  or  S-procedure  argument  [23]. 

Another  upper  bound  to  can  be  computed  by  using  the  implicit  ^-framework 

presented  in  Newlin  [25]  and  Tierno  [36].  From  initial  investigation,  the  implicit  bound  is 
more  conservative  than  the  bound  presented  here. 

3.2.2  Properties  of  JXS 

The  LMI  optimization  in  Theorem  10  can  be  solved  by  using  general  purpose  LMI  solvers, 
but  it  can  be  solved  more  efficiently, 

mf  {7:  M* (I  o  D)M  -  72£>  <  0}  =  (55) 

The  reader  is  refer  ed  to  [30]  for  more  details. 

The  LMI  in  equation  44  is  related  to  the  standard  Jl(M)  LMI,  but  it  is  unclear  whether 
the  new  LMI  is  consistent  with  the  standard  LMI.  Corollary  11  and  Theorem  13  are  presented 


to  validate  that  the  two  bounds  are  consistent. 

Corollary  11  JIS(M)  =  Jl,s{TMT~l)  where  I  o  T  =  T  >  0 
Proof.  The  LMI  in  equation  54  becomes 

T~*M*T*(I  o  D)TMT~l  -  <  0  (56) 

M*T*(I  o  D)TM  —  j2T*BT  <  0  (57) 

Let  T  =  diag[t\,  ■  ■  ■ ,  tn\: 

M*(I  o  D  o  [t*t))M  —  7 2(D  o  (t*t))  <  0  Applying  lemma  1  (58) 

Let  D  —  D  o  (t*t): 

M*(I  o  D)M  —  72Z>  <  0  (59) 

D  >  0  which  follows  from  lemma  2  =7  for  a  given  7,  if  the  LMI  in  equation  56  has  a  solution 
then  the  LMI  in  equation  54  has  a  solution. 

=>  <  Jis(TMT~l)  <  JIS(T~1TMT~1T)  =  Jls(M)  (60) 


o 

In  Corollary  11,  T  enters  in  exactly  like  the  D-scales  from  the  standard  fj,  upper  bound.  This 
implies  that  has  exploited  all  the  structure  in  the  uncertainty  set  that  ~p(M)  has  and 

that  the  additional  freedom  in  the  /JS(M)  LMI  is  related  to  the  structure  of  the  Frobenius 
set  description.  In  the  standard  case,  the  hypercube  does  not  provide  additional  structure  to 
exploit. 

The  following  lemma  is  a  useful  tool  in  relating  Jis(M)  and  /J(M). 
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Lemma  12  D  <  nl  0  D  for  D  >  0 

Proof. 

a(tn)  =  n  and  a(nl )  —  n 

(61) 

=>  nl  —  Jn  >  0 

(62) 

(nl  0  D)  —  D  —  (nl  —  ln)  0  D 

(63) 

Applying  Lemma  2, 

(nl  0  D)  —  D  >0 

(64) 

o 


Theorem  13  <  /x(M)  <  y/nfis(M) 

Proof.  Equations  30  and  54  can  be  rewritten  as 

/Z (M)=  inf  {7:  M*(I  o  D\)M  <  q2(/  o  £>1)}  (65) 

Di>0 

and 

£S(M)  =  inf  {7:  M*(J  o  D2)M  <  72D2}  (66) 

£>2>0 

respectively.  Assume  —  I  °  D\ 

=> Jis(M)<il(M )  (67) 

Apply  Lemma  12: 

M*(I  o  D2)M  <  {Jl s{M))2D2  <  n(7Zs(M))2(J  o  D2)  (68) 

Assume  Di  =  D2' 

=>  <  \fnjis  (M)  (69) 

o 

Theorems  7  and  13  indicate  further  consistency  between  the  spherical  \x  upper  bound,  the 
standard  \x  upper  bound,  and  the  set  containment  shown  in  figure  21.  If  the  relationship 
between  /jls{M)  and  /i(M)  was  not  the  same  as  the  relationship  between  their  respective 
upper  bounds,  then  their  would  be  instances  where  one  of  the  upper  bounds  could  be  trivially 
improved.  For  example,  assume  that  the  relationship  from  Theorem  13  was  replaced  by 
~fis(M)  <  ~p{M)  <  ri]Is(M).  If  both  of  these  upper  bounds  were  computed  with  Jis(M)  =  1 
and  JX(M)  =  n.  Then  latter  upper  bound  would  be  useless,  because  the  former  implies  that 
Ms(Af)  <  1.  This  in  turn  implies  that  //(M)  <  y/n  from  Theorem  7.  So  Jl(M)  could  be 
improved. 

fi(M)  is  equal  to  the  upper  bound  in  the  rank  one  case.  The  exactness  of  the  upper  bound 
for  the  rank  one  case  is  a  direct  consequence  of  convexity.  In  fact,  in  this  case  the  determinant 
is  linear  in  the  uncertain  parameters.  The  same  properties  hold  for  fJ,s(M)  in  the  rank  one 
case. 
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Figure  21:  Cube,  Inscribed,  and  Superscribed  Spheres 


Theorem  14  fis(M )  =  /za(M)  for  rank(M)  =  1 

Proof.  Consider  the  rank  one  matrix  M.  Applying  the  explicit  solution  given  in  equation  55, 
the  optimal  value  of  7  is  7  =  1  lmw|2-  An  explicit  perturbation  that  makes  (I  —  M A) 

singular  is  A  =  (I  o  M')/ y2.  Since  the  Frobenius  norm  of  A  is  I/7  ,  then  7  is  also  a  lower 
bound  for  fis.  Therefore,  fis  =  ~p^  in  the  rank  one  case, 
o 


The  worst  ratio  between  /is(M)  and  /is(M)  is  at  least  yfn.  This  ratio  is  achieved  by 


M  — 


0  1 
I  0  5 


(70) 


for  which  ^  and  fis(M)  —  1. 


3.2.3  Generalizations  of  fis 

The  solution  for  the  upper  bound  of  spherical  [i  can  be  generalized  from  an  uncertainty 
described  by  a  spheres  to  hyperellipsoids.  The  upper  bound  can  also  be  generalized  to  other 
traditional  uncertainty  descriptions  like  repeated  scalars,  real  parameters,  and  full  blocks. 
These  generalizations  can  be  mixed  and  matched  to  form  a  generalization  of  the  LMI  in 
equation  30.  Only  the  resulting  upper  bound  LMIs  will  be  presented.  It  is  relatively  straight 
forward  to  extend  the  derivation  of  the  LMI  in  Theorem  10  to  the  results  presented. 

A  natural  generalization  of  the  hypersphere  is  the  hyperellipsoids, 

8*P5<1,P>0.  (71) 

\x  for  a  hyperellipsoid  is  defined  by 

^  min{(^*P5  :  A  £  As,  det(7  —  MA)  =  0} 
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It  is  trivial  to  generalize  JIS,  equation  44,  when  P  is  a  diagonal  matrix  to  compute  an 
upper  bound ,/Ze(M),  to  elliptical  /r,  fie(M).  For  this  case,  the  principle  axes  of  the  ellipsoid 
are  aligned  with  the  coordinate  axes  in  As.  The  matrix  M  can  be  rescaled  to  create  a  new 
and  equivalent  problem  where  the  uncertainty  can  be  described  by  a  spherical  constraint. 
This  rescaling  is  analogous  to  scaling  a  rectangle  into  a  square.  The  LMI  in  equation  44  can 
be  used  to  compute  an  upper  bound  to  the  original  elliptical  p  problem. 

For  the  general  hyperellipsoid  in  equation  71,  the  proof  in  Theorem  10  can  be  generalized 
to  construct  an  upper  bound  LMI  solution.  The  LMI  for  p  on  the  general  hyperellipsoid  is 
given  by 

pe(M)  <  JIe(M)  =  mf  {7:  M*(P~l  o  D)M  <  72D}.  (73) 

For  Jls(M),  P  —  P~l  =  I  and  equation  73  reduces  to  equation  66.  For  the  LMI  in  equation 
73,  a  closed  form  expression  for  the  optimal  gamma  exists, 

/ Ue{M )  =  yj p((MT  ®  M*)diag(vec{P~1))),  (74) 

the  reader  is  referred  to  [30]. 

Within  the  \x  framework  an  allowable  uncertainty  structure  is  a  repeated  scalar  [], 

^ms:  —  {diag  [8\ICl ,  •  •  • ,  8nICn]  ■  5*  £  C}.  (75) 

This  uncertainty  structure  can  be  fit  into  a  spherical  context.  For  notational  convenience,  a 
matrix 

Q  diag(tCl,-  ■  • ,  lCn),  (76) 

is  introduced.  The  resulting  upper  bound  LMI  is  given  by 

^,Ams(M)  <  7L,  Ams(M)  =  mf{7:  M*(Q  o  D)M  <  72D}.  (77) 

In  equation  77,  Q  is  not  an  invertible  matrix.  Q  enters  the  LMI  as  a  term  of  a  Hadamard 
product,  as  a  result  the  invertibility  of  Q  is  irrelevant.  The  closed  form  solution  is  identical 
to  equation  74  with  P "1  replaced  by  Q. 

The  uncertainty  structure  for  real  parameters  is  described  by, 

Ars:  =  {diag  [Su  •  ■  ■ ,  Sn] :  5{  E  R}.  (78) 

The  additional  freedom  in  the  LMI  in  equation  30  of  the  G-scales  corresponds  to  the  constraint 
that  Si  is  real  and  that  the  input  and  output  signals  of  the  block  must  be  aligned  in  C.  The  G- 
scales  are  independent  of  the  magnitude  and  linkage  between  separate  blocks.  The  resulting 
upper  bound  LMI  is 

Us  Ars  (M)  <  Vs  Ars  (M)  = 

MD>o,G=G*{r-M*(I  O  D)M  +  j((I  O  G)M  -  M* (I  O  G ))  <  72D}.  {  l 

The  uncertainty  structure  for  only  full  block  uncertainty  is 
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Computing  an  upper  bound  for  mu  with  spherical  constraints  a(A^)  which  includes  non-scalar 

Ai, 

X>(Ai)|2<l,  (81) 

i- 1 

is  not  a  natural  extension  of  the  LMI  for  diagonal  uncertainty.  For  diagonal  uncertainty  in 
the  standard  ji  problem  the  structure  of  the  D-scales  can  be  specified  by  Q  o  D  where  D 
is  full.  For  full  blocks  the  resulting  structure  in  the  D-scales  can  not  be  specified  using  the 
Hadamard  product  of  a  constant  matrix  and  a  full  D. 

One  method  to  construct  an  upper  bound  is  to  convert  the  original  spherical  /is,Afs(Af) 
problem  with  full  blocks  into  a  new  problem,  where  M  is 

constructed  from  M.  M  must  be  partitioned  into  M{j  which  is  compatible  with  the  block 
structure  of  the  uncertainty.  What  is  meant  by  this  is  that  M{j  maps  the  output  of  the  jth 
block  of  A  to  the  input  of  the  ith  block  of  A.  M  is  constructed  as  follows, 

Mij  =  a(Mij)  (82) 


The  new  problem  has  a  smaller  dimension  matrix  than  the  original  problem  and  the  upper 
bound  bound  has  a  closed  form  expression,  but  this  bound  will  tend  to  be  conservative  because 
fJ's, Afs(M)  < 

Another  method  is  to  construct  a  larger  problem,  larger  dimension  and  more  blocks  in 
A,  for  which  ^5)Afs(Af)  =  ^,AnewW  and  7*g, Anew(M)  can  be  computed.  The  original 


Figure  22:  Expanded  System 


system  shown  in  figure  20 

has  a  spherical  constraint  on  A.  This  is 

equivalent  to  the  problem 

shown  in  figure  22,  where 

A  =  dio,g(5i Ij\fi ,  •  •  • ,  6^1  €  C 

(83) 

with  the  constraint  that 

\Si\  <  1  ,Vt 

(84) 

and 

A2  =  Aft. 

(85) 
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The  spherical  constraint  can  be  moved  from  A2  to  A1  to  form  an  equivalent  system.  So 
equation  84  is  replaced  by 

n 

J2\Si\2<l  (86) 

i— 1 

and  the  structure  of  A  is  unchanged  but  equation  81  is  replaced  by 

H A?)  <  IV*.  (87) 


The  system  shown  in  figure  22  can  be  put  into  the  original  form  as  shown  in  figure  20.  Let 


By  construction, 


A  ~ 


A1  0 
0  A2 


:  A1  e  A1,  A2  e  A: 


I 


M  = 


0  I 
M  0 


y/vg,  a(M)  =  ^,Afs(M). 


An  upper  bound  of  fj,g  ^  (M)  is  given  by 

inf  <f 

£>i>o,£»2eD2  ^ 


T- 


Q°Di<  i2D2  1 
M*D2M  <j2Dx  J 


(88) 

(89) 

(90) 


(91) 


where  Q  is  similar  to  equation  76,  and  D2  is  the  standard  //  D-scales  for  full  block,  ie  it 
consists  of  Silni  blocks. 
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